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ABSTRACT 


An  automated  general  purpose  system  for  analysis  is 
presented.  This  system,  identified  by  the  acronym,  "MAGIC  III" 
for  Matrix  Analysis  via  Generative  and  Interpretive  Computations, 
is  an  extension  of  the  structural  analysis  capability  available 
in  the  initial  MAGIC  System.  MAGIC  III  provides  a  powerful  frame¬ 
work  for  implementation  of  the  finite  element  analysis  technology 
and  provides  diversified  capability  for  displacement,  stress, 
vibration  and  stability  analyses. 

Additional  elements  have  been  added  to  the  MAGIC  element 
library  in  this  phase  of  MAGIC  development.  These  are  the  solid 
elements;  rectangular  prism,  tetrahedron,  triangular  prism,  sym¬ 
metric  triangular  prism,  and  triangular  ring  (asymmetrical  load¬ 
ing).  Also  included  are  the  symmetric  shear  web  element  and  a 
revised  quadrilateral  thin  shell  element.  The  finite  elements 
listed  include  matrices  for  stiffness,  mass,  prestrain  load, 
thermal  load,  distributed  mechanical  load,pressure  and  stress. 

The  MAGIC  III  System  for  structural  analysis  is  presented 
as  an  integral  part  of  the  overall  design  cycle.  Considerations 
in  this  regard  include,  among  other  things,  preprinted  input  data 
forms,  automated  data  generation,  data  confirmation  features, 
restart  options,  automated  output  data  reduction  and  readable 
output  displays. 

Documentation  of  the  MAGIC  III  System  is  presented  in 
three  parts;  namely,  Volume  I:  Engineer's  Manual,  Volume  II: 

User's  Manual  and  Volume  III:  Programmer's  Manual.  The  subject 
document.  Volume  I  (Engineer's  Manual)  is  an  extension  of  the 
primary  Technical  documents.  Included  are  the  theoretical  develop¬ 
ments  for  the  additional  finite  elements  Included  in  the  MAGIC  III 
System  as  well  as  a  discussion  of  newly  added  computational 
procedures . 
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SECTION  I 


INTRODUCTION 

The  MAGIC  III  System  for  structural  analysis  is  an  extension 
of  the  MAGIC  I  and  MAGIC  II  Systems  reported  in  References  1  to  6. 
All  capabilities  available  in  the  original  systems  have  been 
retained  and  improved  upon.  Extension  of  the  MAGIC  System  has 
been  in  the  following  areas: 

(a)  Incorporation  of  four  (4)  solid  elements 

(1)  Rectangular  Prism 

(2)  Tetrahedron 

(3)  Triangular  Prism 

(4)  Symmetric  Triangular  Prism 

(b)  Incorporation  of  a  Triangular  Cross-Section  Ring  which 
accommodates  asymmetric  mechanical  and  thermal  loading. 

(c)  Incorporation  of  the  Symmetric  Shear  Web  element. 

(d)  Incorporation  of  a  revised  Quadrilateral  Thin  Shell 
element  which  reflects  high  aspect  ratio  usage. 

(e)  Incorporation  of  new  equation  solvers  into  the 
MAGIC  III  System. 

(f)  Inclusion  of  additional  computational  procedures  to 
support  the  analysis  process. 

The  work  reported  herein  is  a  discussion  of  the  extensions 
listed  above.  The  discussion  encompasses  three  volumes  of  which 
this  is  the  first.  This  Volume,  Engineer's  Manual,  (Volume  I)  i;- 
an  addendum  to  the  technical  reports  given  in  References  1  and  4 
and  as  3uch  should  be  used  in  conjunction  with  these  references  to 
effectively  utilize  the  MAGIC  III  System.  The  seconl  Volume,  User's 
Manual,  Reference  7,  includes  detailed  specifications  for  the 
preparation  of  input  data  for  the  additional  elements  in<-iudei  m 
this  third  version  of  MAGIC.  The  last  volume,  Volume  III, 

4 


Programmer's  Manual,  Reference  8,  presents  information  on  the 
organization  of  the  MAGIC  III  System  as  well  as  its  operational 
characteristics . 

Section  II  of  this  report  presents  the  theoretical  basis 
of  the  additional  finite  elements  and  gives  explicit  expressions 
for  their  characteristic  matrices.  These  elements  are: 

(a)  Rectangular  Prism 

(b)  Tetrahedron 

(c)  Triangular  Prism 

(d)  Symmetric  Triangular  Prism 

(e)  Triangular  Cross-Section  Ring  (Asymmetric  Loading) 

(f)  Symmetric  Shear  Web 

(g)  Revised  Quadrilateral  Thin  Shell 

Figures  I-l  to  1-3  depict  these  newly  added  elements  an  well 
as  previously  existing  elements  of  the  MAGIC  System. 

A  discussion  of  new  computational  features  incorporated  into 
the  MAGIC  HI  System  is  given  in  Section  III.  Included  are  a 
discussion  of  the  ANALIC  (Analysis  In  Core)  Module  and  the  out-of¬ 
core  variable  bandwidth  equation  solver  based  on  Cholesky 
triangularization. 

The  body  of  the  technical  report  is  concluded  with  a 
general  retrospective  discussion  in  Section  IV.  An  overview  of  the 
MAGIC  III  System  is  presented. 
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FIGURE  1-1  MAGIC  III  ADDITIONAL  FINITE  ELEMENTS 


FIGURE  1-2  MAGIC  II  ADDITIONAL  FINITE  ELEMENTS 


FIGURE  1-3  MAGIC  FINITE  ELEMENTS 


SECTION  II 


si 

ADDITIONAL  FINITE  ELEMENTS 


A.  INTRODUCTION 

The  MAGIC  III  System  incorporates  seventeen  finite 
elements.  Ten  of  these  elements;  namely  frame,  shear  panel, 
triangular  cross-section  ring,  toroidal  thin  shell  ring, 
quadrilateral  thin  shell,  triangular  thin  shell,  trapezoidal 
cross-section  ring,  quadrilateral  plate,  triangular  plate  and 
incremental  frame  were  available  in  the  initial  MAGIC  and 
MAGIC  II  Systems  and  are  described  in  detail  in  References  1 
and  4. 

5  Seven  additional  elements;  namely  rectangular  prism, 

j  tetrahedron,  triangular  prism,  symmetric  triangular  prism, 

i  symmetric  shear  web,  triangular  cross-section  ring  and  a  revised 

!  quadrilateral  thin  shell  element  have  been  incorporated  into 

*  the  MAGIC  III  System.  Characteristic  matrices  hove  been  derived 

|  for  these  elements  and  include  stiffness,  stress,  prestrain  load, 

■  pressure  load,  thermal  load,  and  consistent  mass  matrices.  The 

j  derivation  of  these  matrices  for  each  finite  element  is  presented 

I  in  the  following  sections. 

i 

I 

1  B.  RECTANGULAR  PRISM  ELEMENT 

I.  Introduction 

The  formulation  of  an  element  stiffness  matrix 
for  the  rectangular  prism  discrete  element  was  first  documented 
in  Reference  9,  and  the  approach  used  here  is  one  of  three 
suggested  therein. 

The  rectangular  prism  element  is  a  powerful  tool 
for  the  analysis  of  solid  structures,  thick  plates,  and  beams. 

It  can  be  used  in  conjunction  with  the  triangular  prism  and 
tetrahedral  discrete  elements  for  the  analysis  of  arbitrary  solid 
geometries,  or  with  plate  elements  for  the  analysis  of  built-up 
regions . 
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An  appropriate  mathematical  model  for  the 
rectangular  prism  discrete  element  is  formulated  on  the  basis 
of  the  variational  principles  of  continuum  mechanics.  From  an 
admissible  assumed  displacement  function  only,  algebraic  expres¬ 
sions  for  various  element  matrices  which  describe  the  mathematical 
behavior  of  the  element  are  derived  by  use  of  the  Lagrange  varia¬ 
tional  equation. 

Consistent  with  the  state  of  the  art,  the  discrete 
element  repre jentation  for  the  subject  element  is  taken  to  consist 
of  algebraic  expressions  for  the  following  matrices : 


a. 

Stiffness 

[X] 

b. 

Stress 

ES] 

c. 

Prestrain  Load 

{Fe} 

d. 

Thermal  Load 

Eft) 

e. 

Consistent  Mass 

EM] 

f. 

Pressure  Load 

EFp} 

These  matrices  arise  as  coefficient  matrices  in 
the  generalized  form  of  the  Lagrange  equations.  The  1'orm  of  these 
equations,  necessary  for  the  complete  element  representation  listed 
above,  are: 


where: 


4*  Vi 

%  generalized  displacement 

Hr 

4*  Vi 

‘  %  generalized  velocity 

Hr 

'v  total  potential  energy 
3^  'v  total  kinetic  energy 
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II.  Geometry 

Figure  II-l  depicts  the  geometry  of  the  rectangular 
prism  element.  Also  shown  are  the  local  and  global  axes  systems; 
namely,  local  x,  y,  z  and  global  X,  Y,  Z.  The  local  axes  are  fixed 
at  the  centroid  of  the  element.  Use  of  vector  analysis  permits 
definition  of  the  dimensions  of  the  prism  to  be: 


a  =  1/2  |rx| 
b  =  1/2  i?y| 
c  =  1/2  |r2| 


(2) 

(3) 

(4) 


Where : 

l?xl  ■ 

IV?8l  = 

[<W 

II 

IV*8l  = 

[(x7-x8) 

\rz\  ‘ 

!r5-fg|  =  | 

>5-x8) 

1/2 

(Y4-Y8>2  + 

(zrz8)2] 

(5) 

1/2 

<Vy  8)2  + 

(z7-z8)2] 

(6) 

0 

1/2 

(7) 

«5-V  + 

(z5-z8)2] 

The  quantities  r^,  r^,  and  rg  are  vectors  emanating 

from  the  origin  of  the  global  axes  to  prism  grid  points 
4,  5,  7  and  8.  The  vectors  r  ,  r  ,  r  form  a  mutually 

*»  J  " 

orthogonal  set  (see  Figure  II-l). 


A  rotational  and  translational  transformation 
matrix  from  local  to  global  coordinates  is  formed  using  these 
vectors.  This  transformation  Is  given  below. 
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f 


(*<*>!  »  [V^X<g))-{Xo)(S>)^ 


where : 


{X(^)}T  _  LxU)j  ,(Oj  are  the  local  coordinates 

{X^)T  =  |x(gJ,  Y(g\  Z(g^J  are  the  global  coordinates 


{Xc(g)}T=  |x£g),  ?cS\  z£g)J  are  the  centroidal  global 


coordinates 


tV  =  pv  %■  e*3 


is  the  matrix  of 
6yl*  Sy2*  Sy3  direction  cosines 


ez  j  ez  *  e* 

zi  za  "3 


-k  (x8  '  xa)*  V 

4(VX8>>  % 


—  (X5~X8^»  ez 
2c  p  0  z2 


■ia  VV  •  V 

-k  <vv>  V3- 

4r  <W>  V  * 


-k  <Z^8> 

IT  (V*8> 


2c  ( Z^-Zg ) 


The  transformation  matrix  [T  ^3  is  Used  not  only  for  coordinate 
transformations  but  deformation  transformations  also. 
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III.  Assumed  Displacement  Functions 

A  structural  element  Is  mathematically  dis¬ 
cretized  into  a  finite  number  of  displacement  degrees  of  free¬ 
dom  by  the  assumption  of  displacement  mode  shapes.  For  the 
simple  geometry  of  the  rectangular  prism  element,  trilinear 
Lagrangian  interpolation  formulas  are  constructed.  The  dis¬ 
placement  is  given  by 


5(^)(x,y,z)  =  l/8abc  )  ~f!f2V3(' * ) 


(d) 


-  f,f. 


-f,fnf06,:'‘J  '  +  f,  f0f-6_ 
X  d  i  D  X  d  5  ( 


WaV3^ 


(9) 


where 


fx  =  (x+a)  f2  =  (y+b)  f^  =  (z+c) 

f-L  =  (x-a)  f2  =  (y-b)  f3  *  (z-c) 


(10) 


and  a,  b,  and  c  are  the  half-dimensions  of  the  prism  as  shown 
in  Figure  II-l. 

Note  that  ^  >  k  *  1>  2,  .  .  .  .  8  are  the 

grid  point  displacements  where  j  =1,  2,  3  corresponds  to  the 
u,  v,  and  w  displacements. 

Equation  (9)  can  be  written  in  matrix  form  as 


6^ )  =  l_  Bj  {5^  ^ )  *  8  abc 


(11) 


where 


(B)  =  L-f1f2f3»  ~^*i^>2^>3 *  ^1^2 ^3 *  ^*1^2 ^3*  "fif2f3< 

?lV3*  “^2^3  J 


(12) 


and 


{5^))T=  bl(J)> 

(13) 
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It  is  instructive  to  examine  the  nature  of  these 
assumed  displacement  functions  by  considering  the  allowable  de¬ 
formation  of  .-ach  face  of  the  prism.  For  example j  the  displace¬ 
ment  of  the  planes  x  =  a,  y  =  b,  and  z  =  c  can  be  written  for 
x  =  a 

S^(a,y,z)  =  (k] yz  +  k^y  +  k^z  +  +  «2^  - 

+  64^)  (1*0 


for  y  =  b 

6(<})(x,b,z)  =  (k2  xz  +  k2  x  +  +  ki|)(62(t})  -  S3(,}) 

66(,3)  +  (li>) 

for  z  =  c 

6^(x,y,c)  =  (k1  xy  +  k2  x  +  k3  y  +  k^-fii^)  +  5^*^ 

+  6  6  <3))  (16) 

5  £ 

and  similarly  for  x  »  -a,  y  =  -b,  z  =  -c.  It  is  noted  that 
the  k1  are  arbitrary  constants.  Referring  to  Figure  II-l, 
it  is  seen  that  the  displacements  on  these  planes  are  functions 
only  of  the  displacements  of  the  gridpoints  defining  the  planes. 
Hence,  the  assumed  functions  are  admissible  in  that  they  satisfy 
the  requirements  of  displacement  continuity  along  interelement 
boundaries.  Due  to  the  assumption  of  linear  interpolation 
formulas,  the  edges  of  the  prism  remain  linear  in  deformation. 

A  direct  consequence  of  the  above  observations  is  that  although 
a  single  element  rnaywarp  under  a  force-couple,  it  may  not  bend 
under  any  conditions. 

The  foregoing  assumed  displacement  functions 
lead  to  three  translational  displacement  degrees  of  freedom 
at  each  of  the  eight  corner  gridpoints;  thus  the  complete 
element  deformation  is  described  by  twenty-four  displacement 
degrees  of  freedom. 


* 


I 
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It  is  instructive  to  examine  the  nature  of  these 
assumed  displacement  functions  by  considering  the  allowable  de¬ 
formation  of  each  face  of  the  prism.  For  example,  the  displace¬ 
ment  of  the  planes  x  =  a,  y  =  b,  and  z  =  c  can  be  written  for 
x  =  a 

fi(J)(a,y,z)  =  (k^yz  +  k^y  +  lc^z  +  +  S2^  -  6^^ 

+  64(J))  (1*0 


for  y  =  b 

6^(x,b,z)  =  (k1  xz  +  k2  x  +  k^2  +  k^)  (62  ^  ^  -  63^ 

_  S6(J)  +  57(J))  (lb) 

for  z  =  c 

6^^(x,y,c)  =  (k^  xy  +  k2  x  +  k3  y  +  kjj)(-6i^)  + 

+  6  6  ) )  (16) 
b  6 

and  similarly  for  x  =  -a,  y  =  -b,  z  =  -c.  It  is  noted  that 
the  k1  are  arbitrary  constants.  Referring  to  Figure  II-l, 
it  is  seen  that  the  displacements  on  these  planes  are  functions 
only  of  the  displacements  of  the  gridpoints  defining  the  planes. 
Hence,  the  assumed  functions  are  admissible  in  that  they  satisfy 
the  requirements  of  displacement  continuity  along  interelement 
boundaries.  Due  to  the  assumption  of  linear  interpolation 
formulas,  the  edges  of  the  prism  remain  linear  in  deformation. 

A  direct  consequence  of  the  above  observations  is  that  although 
a  single  element  maywarp  under  a  force-couple,  it  may  not  bend 
under  any  conditions. 

The  foregoing  assumed  displacement  functions 
lead  to  three  translational  displacement  degrees  of  freedom 
at  each  of  the  eight  corner  gridpoints;  thus  the  complete 
element  deformation  is  described  by  twenty-four  displacement 
degrees  of  freedom. 
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The  definition  of  assumed  displacement  functions 
permits  the  derivation  of  the  strain-displacement  relationships. 
The  element  strain  components  are  expressed  as  functions  of  the 
assumed  displacement  modes  by 


e 

X 

-  i(i>  -  ■ 

y  A 

S  <5^ 

3x 

(17) 

e 

»  fi(2)  =  . 

3S(2) 

(18) 

y 

,y 

3y 

£  : 
z 

-  6(3)  .  . 

jZ 

as  ( 3) 

3z 

(19) 

exy 

=  6(l)  + 

fi(2>  _  36^ , 

36(2) 

(20) 

y  y 

jX  "  3 

9y 

3x 

eyz 

=  +  <5 

>z 

<3>  -  36<2) 

>y  -  36  + 

3z 

36<3> 

3y 

(21) 

ezx 

=  6^  '  + 
>z 

“(1)  + 
3z 

36(3) 

3x 

(22) 

Performing  the  necessary 
displacement  functions  yields: 

differentiations 

«<J> 

1 

lp*J<«£m 

(23) 

8  abc 

,y 

1 

LDyJtsk3>> 

8  abc 

(24) 

6(j ) 
z 

B  1 

ug<4J>> 

8  abc 

(25) 
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where 


'V 


<v 


f2f3*  ^2**3»  ^2f3*  "f2f 3 

f2^3*  -f2f3^ 

L_rif 3 *  flf3J  “f1^3*  fl*T  ^lf3*  ~*lf3* 
*1*3*  -f1f3j 


(2b) 


(27) 


(Dz>  L  ^1^2*  ^1^2*  “^1^*2*  ^1^2*  ^1^2 

-f^J  (28) 


and 


{s£J)}T=  L6iJ),  ^2(J),  63(J),  6„(j),  6S(J), 


«.<*>,  «7(J},  «.(J>j 


(29) 


These  equations  are  assembled  into  a  single 
matrix  expression  relating  the  strain  components  to  the  dis¬ 
placement  degrees  of  freedom. 


where 


<c)  =  -5-V  [D] 

8  abc 


(30) 


eyz*  ezxJ 


(6^>}T 


(3D 


(32) 


Potential  Enerp 


The  potential  energy  of  the  element  is 


6  =  u-W 

P 


where 


lei 

U  =  j  f  [dej  {0} 


v  “'{0} 


W  are  external  work  contributions 


(e)  =  Lex»  ey»  exy>  eyz»  ezxJ 


{o}A  =|o, 0.0,0  ,  O  ,  0  I  (38) 

L  x*  y*  z*  xys  yz’  zxJ  •  J  ' 


Linear  elastic  materials  behaviour  is  assumed  from  an 
initial  state  of  strain  {e}  to  a  final  state  of  stress  {0} 
and  strain  {e}.  From  the  generalized  Hooke's  Law, 


{0}  =  CE][{e}  -  {?}} 


where  [E]  is  the  symmetric  matrix  of  elastic  constants 
which,  for  three-dimensional  orthotropic  material,  can  be 


written 


Ex(1-VV-  E*(V  +  VwV 


[ E ] = 1/ A 


,  E  ( 1— v  v  ),  E  (v  +v  v  )  0  0  0  ( 40) 

*  yv  zx  xz  ’  y  zy  xy  zx  u»  u  v  ' 


,  E  (1-v  v  ) 
*  zN  xy  yx' 


-Symmetrlc- 


0,  0,  0 

.4Gxy>  0 
*aGyz  ■  0 
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where 


U  =  J  (1/2  J_eJ  [E]  {e}  -  IcJ  [E]  {im  .  (*<3) 

V 

This  is  the  desired  form  of  potential 

energy  4>  . 

r 


V. 


Element  Static  Matrices 


5.1  Introduction 

To  effect  the  discretization  of  the 
element  the  assumed  displacement  functions  are  introduced 
into  the  potential  energy  function  which  in  turn  is  substi¬ 
tuted  into  the  Lagrange  equations  to  yield  element  matrices 
with  respect  to  grid  point  displacement  degrees  of  freedom. 

An  exception  is  the  element  stress  matrix  which  is  derived 
from  strain-displacement  and  stress-strain  relationships. 

5.2  Stiffness  Matrix 

The  energy  contribution  to  the  linear 
elastic  stiffness  is  given  by 

*  *  1/2  f  [ej  [E]  {e}  dV  .  (44) 

v 

Recalling  the  strain-displacement  relations 

{e}  ,  _L_  [D]  {6(J)>  (30) 

8  abc 

and  substituting  these  into  Equation  (44)  yields 

%  ■  I  <BS5J>2  JL«U)JMT  CEJ  tD]  <«(J))dV  .  (45) 

V 


Performing  the  matrix  multiplication,  and  noting  that  the 
grid  point  displacement  degrees  of  freedom  are  not  functions 
of  the  local  coordinate  system,  the  potential  energy  function 
may  be  expressed  as  shown  by  Equation  46.  Taking  the  first 
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AP  rV|{xa> 


variation  of  the  potential  energy  function  with  respect  to 
displacement  degrees  of  freedom  (as  showr.  by  the  first  term 
of  Equation  (1)),  yields  the  element  stiffness  matrix  [K] 
referenced  to  local  grid  point  displacements.  This  matrix 
is  depicted  in  Equation  (47). 


The  matrix  products  appearing  as  integrands  in 
Equation  (48)  lead  to  integrations  of  the  following  general 
form, 


-  ■/// 
-a  -b  -< 


(x+a)m  (x-a)n  (y+b)p  (y-b)q  (z+c)r  (z-c)sdzdydx 

(49) 


Since  the  limits  of  integration  arc  constants.  Equation 
(49)  can  be  written  equivalently  as 


(x+a)m(x-a)n 


(y+b)p  (y-b)q  dy 


c 

(50) 

-c 

These  definite  integrals  are  readily  evaluated  by  integra¬ 
tion  by  parts  and  the  [I^j]  matrices  are  expressed  in 
Equations  (51)  -  (56). 


(z+c)r  (z-c)s  dz 
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i£  a2b2c3 
3 


-2  2  1 
-2  2  1 
-112 
-112 
2  -2  -1 
2  -2  -1 
1  -1  -2 
1  -1  -2 


-1 


-1 

-2 

-2 

1 

1 

2 

2 


-2  2  1-1 

-2  2  1-1 

-112-2 
-112-2 
2-2-1  1 

2-2-1  1 

1-1-2  2 

1-1-2  2 


CV 


(54) 
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3 


aW 


-2-22  2 

2  2-2  -2 

2  2-2  -2 

-2-22  2 

-1-11  1 

1  1-1  -1 

1  1-1  -1 

-1-11  1 


-1-11  1 

1  1-1  -1 

1  1-1  -1 

-1-11  1 

-2-22  2 

2  2-2  -2 

2  2-2  -2 

-2-22  2 


[Iyz] 


(55) 

'! 


16  a2b3c2 
3 


2  11  2 

12  2  1 

-1  -2  -2  -1 

-2  -1  -1  -2 

2  11  2 

12  2  1 

-1  -2  -2  -1 

-2  -1  -1  -2 


-2  -1  -1  -2 

-1  -2  -2  -1 

12  2  1 

2  11  2 

-2  -1  -1  -2 

-1  -2  -2  -1 

12  2  1 

2  17  2 


(56) 

i  1 

i 
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The  potential  energy  function  given  by 
Equation  (46)  is  referenced  to  local  gridpoint  displacements. 
These  displacements  must  first  be  reordered  to  be  compatible 
with  the  MAGIC  III  ordering  system  and  then  be  transformed  to 
global  displacements.  The  former  is  accomplished  through  use  of 
the  transformation  given  below: 


{S(J)}  =  [T]  {6(J)}  (57) 


The  transformation  from  local  to  global 
deformation  is  derived  through  the  use  of  Equation  (8),  thus 


{ 


{6(J)>  =  [T  - ]  {A} 


wnere 


{A}T  »  JUi.  Vw  Wx.  U2-  V2^  W2^  ....  U„V„KtJ 


The  U,  V,  W 
Also, 


deformations  are  defined  in  Figure  II-l. 


=  [T*  3 


% 

['T  ] 

°9. 


CV 


W^3 


CV 


^  -i 


tv3 


The  [T  .]  matrix  is  defined  by  Equation  (8). 

gx, 


/ 
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Use  of  Equations  (57)  and  (58)  in  the 
potential  energy  equation.  Equation  (46),  yields 

*p  =  1  (— ]T  [T]T  [K]  [T][T  ]  {  A) .  (59) 

r  ~T  8  abc  E*  g* 


Taking  the  first  variation  of  the  potential  energy  with 
respect  to  the  displacement  degrees  of  freedom  {A}  (as  shown 
by  the  first  term  in  Equation  (1)) yields  the  element  stiffness 
matrix  [K] referenced  to  global  gridpoint  displacements.  This 
matrix  is  depicted  by  Equation  (60). 

[K]  =  (— —  )2  [T.]T  [T]T  [K]  [T]  [T  ]  .  (60) 

8  abc  g* 

5.3  Stress  Matrices 

The  element  stress  matrix  follows  as  a 
direct  consequence  of  the  strain-displacement  and  stress-strain 
relationships.  Recalling  that 

{0}  =  CE]{(e}  -  .{e}}  (39) 

where  {£}  is  a  column  of  either  mechanical  prestrain  or 
thermal  prestrain  or  both.  Also  recalling  that 

{e}  =  (__!_)  [D]  (30) 

8  abc 

it  follows  that 

(a)  =  CE]  CD]  {6(J)}  -  [E]  {£}  -  [E]{ea>  •  (6l) 

Use  of  Equations  (57)  and  (58)  in  Equation  (6l)  yields: 

{0}  =  (o~~)  CE]  CD]  [T]  [T  ]  (A)  -  [E]  {e}  -CE]{f>(62) 

8  abc  g* 
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The  distribution  of  prestrain  throughout 
the  element  is  assumed  to  be  of  the  same  functional  form  as  the 
displacement  mode  shapes;  i.e.,  an  interpolation  between  grid 
point  prestrain  components.  Thus, 

eU>  =  (  -1. .  -  )  [Bj  {F(j)}jk  =  1,  2,  8  (63) 

8  abc 

with  j  =1,  2,  .  6  corresponding  to  e  ,  e  ,  e  ,  e  , 

x  y  z  xy 

e  •  i 
yz’  zx  . 

The  vector  {e}  now  becomes 

{e>  =  (q-~-  )  [B]  {e(J)}  (6U) 


where 
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The  vector  {e0},  the  prestrain  due  to  thermal  effects,  can 
be  written  as 


where 


{£*}  =  (_J_)  [B]  [a]  {AT.  }  (65) 

8  abc  K 


ax 

[i] 

ay 

[i] 

az 

[i] 

a 


x> 


a. 


[I]  is  an  (8  x  8)  identity  matrix  and 
are  coefficients  of  thermal  expansion. 


(ATk}T  =  LTk  -  Tj,  k  =  1,  2,  3  8 


x.i_ 

Tk  is  the  temperature  at  the  ic^  grid  point  and  Tq  is  the 
element  reference  temperature. 

Equation  (62)  can  now  be  rewritten  as 


follows: 


{a}  =  [S]  {A}  -  (s)  -  (s‘} 


(66) 


Each  term  in  Equation  (56)  has  a  particular  significance. 

The  matrix 

tS]  =  (FlbI}  CE]  CD]  [T]  CTg£3  (67) 

yields  the  element  apparent  stresses  due  to  displacements 
of  the  element  and  is  referred  to  as  the  element  stress  matrix. 

The  matrix 

(s)  =  (~— )  [E]  [B]  {£  (J)}  (68) 

o  abc  K 


28 


yields  stresses  due  to  the  prestrain  state  within  the  element 
and  is  referred  to  as  the  element  stress  matrix  due  to  pre¬ 
strain.  The  matrix 

)  *  <B“ibc  }  CE]  Cb]  [“3  UV  (69) 

yields  stresses  due  to  a  temperature  gradient  within  the 
element  and  is  referred  to  as  the  element  thermal  stress 
matrix. 

It  is  noted  that  the  assumption  made  in 
Equation  (63)  is  invalid  for  a  constant  temperature  and  prer, train 
distribution  throughout  the  element  since  this  assumption  pro¬ 
duces  zero  prestrain  and  zero  thermal  forces.  Thus  for  a 
constant  prestrain  and  temperature  distribution,  the  following 
equations  replace  Equations  (68)  and  (69). 
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Performing  the  multiplication  and  integration  gives 


=  l6(J)J  tPg} 


where 


{ F- }  : 
e 

(8  abc} 

C?3  (eK 

(J)} 

”En 

CIXB3, 

E12 

[Ixb3> 

E13 

C 1  yb  3 » 

euh 

[I 

yB3,  °» 

e66  [iZB 

LP3- 

E12 

C*YB3  * 

E22 

[lyB3> 

E23 

C I Yg] > 

Elfl 1 

Cl 

XB  3  ’  E 

55  [IZB3, 

E13 

C^  ZB3 » 

E23 

^ZB3, 

E33 

tJZB3 » 

°.  E55 

[Iyb3> 

E66  CIXB 

^XB3  =  1 

{DX 

}  LbJ  dV  »  a  [ixx] 

'V 


[IyB3  =  [  {Dy>  LbJ  dv  =  b  [Iyy] 

CiZb3  =  f  <Dz}  lAl  dV  51  c  ^zz3 
Jv 

Transformation  of  |_6  to  global  gridpoint  deformations 
through  use  of  Equations  (57)  and  (58)  and  differentiation 
of  the  result  with  respect  to  the  gridpoint  deformations 
yields  the  prestrain  load  vector 

(Pg)  =  [TgA]T  [T]T  {F- }. 
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(69) 

] 

0  ( 69  A ) 

] 


(70) 


For  a  constant  prestrain  throughout  the  element 


*-  =  /  [6(J)JC/D]T  [E]  {e}  dV  (71; 

e  8  abc  YJ 

Thus : 

{F,}  =  CP]  { e }  (72) 

where  [P]  is  given  by  Equation  (69B)  and  the  load  vector 
is  given  by  Equation  (70). 

5.5  Thermal  Load  Matrix 

The  thermal  load  matrix  is  a  special  case 
of  the  prestrain  load  matrix.  Substitution  of  the  thermal 
strain 


=  Ca]{ATk) 


into  Equation  (69)  defines  a  load  matrix 

‘V  -  <rk)2  [p,]  *  {4V 


where 


CP  ]  = 


E11  ^XB-^  E12  ^XB^  *  E13  ^XB3 


E12  ^YB^ 3  E22  ^Yll*  E23  ^YB^ 
E13  ^ZB^*  E23  ^  IZB^  *  E33  ^ZB^ 


(73) 


(7^) 


31 


The  final  thermal  load  matrix  is  given  by 


{Fa}  =  [Tg,]T  [T]T  {Pa} 


(75) 


For  a  constant  temperature  distribution 
throughout  the  element,  the  thermal  strains  are  given  by 


{e}“  =  ATave>  («} 


(76) 


where 


AT  =  T  -  T 
fliave.  ave.  o 


T  _  =  T,  +  T0  +  ...  +  T« 

ave .12  o 


Thus,  substitution  of  Equation  (76)  into  Equation  (71) 
defines  a  load  matrix 


<*«>  -  “ave.  ["P'3  <«» 


(77) 


1/ 

The  [P1]  matrix  is  obtained  by  taking  only  the  first  three 

A# 

columns  of  the  [P]  matrix. 
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5.6  Pressure  Load  Matrix 

'..‘he  pressure  load  matrix  is  derived  on 

the  basis  of  constant  pressure  on  each  face  of  the  rectangular 

prism  element.  Thus,  the  total  work  W  due  to  the  pressure 

P 

loads  is  the  sum  of  the  work  done  on  each  face. 


Wp  =  W1234  +  W5678  +  •  •  •  +  woii 


3478 


The  subscripts  denote  a  face  of  the  prism  (see  Figure  II-l). 
Now 


“lsst  J 


■>1234  6<1>  dA 


Recalling  Equation  (11),  Equation  (79)  can  be  written  as 


a  b 


'1234 


_1 _  |bI  {6 ^ 1  ^ }  dxdy 


8  abc 


-a  -b 


Performing  the  indicated  integration  yields 

W1234  =  p1234  bc  L5^3)’  62(1)’  63  (1)>  6,(1)J{1) 

Additional  integrations  of  the  form  shown  by  Equations  (79) 
and  (80)  for  the  remaining  faces  yields 


(3>,  «.C3>, 

67C3),  6»C3)J{1) 

(82) 

(2)  ,  (2) 
y  °  3  y 

6e<2>,  67<2U^> 

(83) 

<2\  6„(2), 

6s(2),  68(2)J{1) 

(84) 
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^1256  =  ^1256  L6»(  ^  *  fi  s  ^  fi«.^  (1) 


(85) 


(i)  (»)  (O  0) 

^5678  =  ^5678  !_?*  >  fie  J  (1)  (86) 


Equations  (81)  to  (86)  are  combined  to  yield  total  /fork 

w  -  L«eU)J  t^p!  -  UJ  [Tg,]T  [T]T  IPp) 


(87) 


through  use  of  Equations  (57)  and  (58)  and 


'v  T  . 

{ Pp )  -  [_F  1  F 1  F 1  Pi  F  2  Pa  P2  P2  F3  F<*  Ph  F  3 

*  ))>>))))})>) 


P3,  F«.,  Fh,  F,,  FS,  fs,  F6,  F6,  Pi,  Ps,  P6,  FeJ 


Where: 


F 1  =  P1231  be,  F2  *  p 5 s 7 1  be,  Ps  -  Pit58  ac 
FH  *  p2367  ac,  F5  -  pl  2  56  ab,  P6  *  p3l(78  ab  . 


Taking  the  first  variation  of  Wp  in  accordance  with 
Equation  (1)  yields  the  pressure  load  matrix 


(Fp)  =  [Tgi]T  [T]T  {Pp} 


(88) 
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VI '  Kinetic  Energy  -  Mass  Matrix 

system  assumin JT  “T*”  ^  8  disorete  “^anlcal 

*  ’  assuming  a  constant  mass  density  p,  can  be  written 

=  o/2  J  IqJ  [J]  (q)  av  (89j 


♦k  - 


veiorciUeVlichldrtity  matrlX  Wher6  W  are  Valued 

locities  which,  from  the  assumed  displacement  modes,  are: 


{iiT  -(!)  -(3) 

iqi  "  L6  ,6,6  J 


Substituting  Equation  (90)  into  Equation  (89) 


(90) 


yields 


$K  =  p/2  J  \}{  \  6(2),  6(3)J  ClJ 

V 


«(l) 

6^ 

*(3> 


dV 


(91) 


or 


*K  -  P/2  /  [(6(,))2  +  (6(2!)2  t  ({<3^3  dv 
V 

Recalling  the  displacement  mode  shapes 

<(,)  -  rV-  LSJ  u<J>> 

o  abc  k 

and  differentiating  gives 

«0)  -  tV  w  &*>} 

0  abc  k 


(92) 


(11) 


(93) 
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Thence, 

(gU))2  =  (__L_  )2  |i(j)J  Cm]  {6v(j)>  <9*0 

8  abc 

where 

'V  %  i 

[m]  =  {B)  [Bj  • 

Substituting  Equation  (9*0  into  Equation  (92),  for 
j  =  1,  2,  3,  yields 


*K 


,  2  r  i\  .(2)  .  (3), 

=  (_~1 — )  p/2  j  [S(  1| 

8  abc  J 


V 


«0>" 

Cm] 

Cm] 

■(2) 
6[  } 

;(3) 

Cm] 

6 

i 

dV 


(95) 


or 

<j>  -  [j(j)J  CM3  {6(J)} 

K  2 


(96) 


where 


CM] 


8  abc 


>2  P 


J 

v 


f\, 

[m] 


(97) 
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Performing,  as  previously,  the  integration  indicated  in 
Equation  (97)  gives 


«  “  f  <rk)7  ["]  dv  = 


,8  abc 
27  p 


1/2  1 
1/4  1/2 


SYMM. 


1/2  1/4  1/2  1 

1/2  1/4  1/8  1/4  1 

1/4  1/2  1/4  1/8  1/2  1 

1/8  1/4  i/2  1/4  1/4  1/2 

1/4  1/8  1/4  1/2  1/2  l/u 


1/2  1 


hence 

[M]  *  [m],  0,  0 

0,  [m],  0  (99 

0,  0,  [m] 

Reordering  and  transform-'  ;o  global  deformations  through  use 
of  Equations  (57)  and  (d^;  permits  the  kinetic  energy  to  be 
written  as: 


4>  =  -lii~  CT  .1T  [T]T  [ft]  [T]  [T  J  U) 


(100) 


Taking  the  first  variation  of  with  respect  to  velocities 
and  differentiating  once  with  respect  to  time,  as  shown  in 
Equation  (1),  yields  the  desired  consistent  mass  matrix 
referenced  to  global  gridpoint  displacements.  Thus, 

CM]  =  [TgA]T  [T]T  [M]  [T]  [Tga]  .  < 


(101) 
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TETRAHEDRON  ELEMENT 


C. 

I.  Introduction 

The  stiffness  matrix  for  the  tetrahedron  element  was 
first  derived  and  presented  in  References  10  and  11  respectively. 
Later  these  relationships  were  reviewed  and  a  consistent 
mass  matrix  was  reported  in  Reference  12.  These  formulations 
have  been  extended  in  MAGIC  III  to  include  stress  fsj  , 
prestrain  load  £f}  ,  thermal  load  ^F^  ,  and  pressure  load  ^Fp] 

matrices.  These  matrices  were  formulated  on  the  basis  of  the 
variational  principles  of  continuum  mechanics.  The  material  that 
follows  summarizes  the  derivation  of  all  the  element  matrices 
mentioned  above. 

A  linear  polynomial  is  assumed  for  each  of  the 
three  displacement  modes.  These  mode  shapes  lead  to  a  total  of 
twelve  undetermined  coefficients  for  the  element  which  are  chosen 
to  correspond  to  three  translational  displacement  degrees  of 
freedom  at  each  of  the  four  vertices  of  the  element.  The  nature 
of  the  assumed  displacement  modes  is  such  that  the  strains 
throughout  the  element  are  constant. 

The  tetrahedron  element  can  be  used  to  analyze 
solid  structures  such  as  thick  plates  and  beams.  It  can  also  be 
used  in  conjunction  with  the  rectangular  prism  and  triangular 
prism  solid  elements  and  in  fact  is  used  to  generate  the 
triangular  prism  element. 

II.  Geometry 

Figure  II-2  depicts  the  geometry  of  the  tetrahedron 
element.  The  local  axes  system  x,  y,  z,  and  global  system  X,  Y,  Z 
are  also  shown.  The  local  axes  are  fixed  at  element  gridpoint  one 
with  the  positive  x  axis  directed  along  side  one-three  as  shown. 
The  global  coordinates  of  each  grid  point  are  given  as 
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AP  - 

*432 

3ViBi,ii 

.xz 

,A432  “ 

3vlEt,ll 

aP  = 
,*432 

3V|B6>1( 

(6) 

aP  = 
431 

3V|Bij2| 

aP  = 

,*431 

3V|b„i2| 

axy  - 
,*431 

3VIB6.2I 

(7) 

aP  = 
*421 

3V|b1j3| 

aP  = 

,  *421 

3V|b4;3I 

axy  - 
A421  ~ 

3V|B6;3I 

(8) 

AYZ  = 
*321 

3VlBj.nl 

AXZ  = 

,  *321 

3v|b4j1|| 

AXY  = 
,*321 

3V|B6j4| 

(9) 
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The  terms  | ^ |  ,  jB^  2|  ,  etc.  represent  the  absolute  value 

of  elements  (1,1)  and  (1,2),  etc.  in  the  LB]  matrix  which  relates 
strains  to  displacements  as  shown  in  Equation  (19)  of  this  Section. 

A  rotational  and  translational  transformation 
matrix  from  global  to  local  coordinates  is  formed  thru  definition 
of  position  vectors  emanating  fro«  the  origin  of  the  global  axes 
systpw  to  element  grid  points  1,  2  ami  3.  This  transformation 
matrix  is 

{*<«}  -  [TgJt3  -  {xi(g)*i  UQ) 


where 

(x<*V  »  [jr^*3,  y^*3"  2^3j  are  the  local  coordinates. 
{x(g),T  .  jU)(  Z(E)J 

are  the  global  coordinates. 

{Xi(«>>T  -  Lxi<g>,  *.ts).  zf*>J  are  the  global  coordinates 

of  gridpoint  one. 


*11 

*12 

e13 

Ti 7F  a 

liTT" 

*  itrr 

e21 

e22 

e?3 

"W  * 

lei  1 

*  lexl 

e31 

e32 

e33 

leal  * 

|e,| 

l®i  1 

t2 


III.  Assumed  Displacement  Functions,  Strain- 
Displacement 

The  assumed  displacement  functions  in  the 

/ 

global  coordinate  system  are 

ux  «  C^X  +  C2  Y  +  C3Z  +  C4 


(11) 


Uy  *  C5X  +  C6Y  +  C7Z  +  C8 


u. 


V  +  °10Y  +  C11Z  +  C12 


(12) 

(13) 


where  Ux,  Uy,  Uz  are  the  deformations  of  the  element 
along  the  global  X,  Y»  and  Z  axes. 

Evaluation  of  Equation  (11)  at  the  four  gridpoints  yields: 


where 


■  (uY }  -  cri  (o 

xi 


{Ux  }  *  LUx  »  Ux  >  Ux  *  Ux,J 

Aj  Jig  A3  A  || 

(c}T  -  lclt  c2,  c3,  C^J 


(14) 


cn 


1  X1  Y1  Z1 
1  X2  y2  Z2 

1  X3  y3  Z3 
L1  X4  y4  y4 


Thus  {C}  -  [P]"1  CUX^}  &  Ux  »  l_l,  X,  Y,  Zj  [ r3“A  {Ux^}  (15) 


1-1 


Likewise,  U  «  jjL ,  X,  Y,  Zj  [H"1  (Uy  } 
y  J  i 


(16) 


uz  -  (i»  x,  y,  zj  [rr1{u  >  •  (17) 

Equations  (15)  to  (17)  are  used  to  derive  element  matrices. 
Note  that  the  displacements  functions  are  written  in  terms 
of  global  coordinates  and  displacements. 

44 


Definition  of  the  assumed  displacement 
functions  permits  derivation  of  the  strain-displacement 
relations.  The  element  strain  components  are 


{e}  = 


~ex  ‘ 

"  9Ux/3x 

ey 

3Uy/3y 

ez 

3U  /3z 

Z 

Yxy 

= 

3U„/3y  +  3U  /3x 

A  y 

Yyz 

3U  /3z  +  3U  /3y 
y  L 

Yxz 

3UZ/3X  +  3UX/3Z 

=  [B]  {U} 


(18) 


where 


(u!  *  L«v  »v  ^  “y3»  «yi|. 

V  V  V  V 

The  [B]  matrix  is  constructed  from  the  [r]-1  matrix  as  follows 

- Zeroes  •  - 


[B]  = 


Row  2  of  Cr]-1, 


-Zeroes 


*;Row  3  of  [rr*- 


—  Zeroes  - 
,-l 


- Zeroes - 

■Row  4  of  [r]-1 


Row  3  of  Cr]"1,  Row  2  of  [r]"1  } Row  3  of  [r]' 
Row  4  of  - Zeroes - »-3Row  2  of  [T] 


-1 


(19) 
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IV 


Potential  Energy 


The  desired  form  of  the  potential  energy  is 

U  -  j"  (|  |ej  [E]  {e}  -  [_ej  CE]  (e))  dV  (20) 

V 

which  was  derived  in  Section  II-B.IV  of  this  report.  The  matrix 
[E]  is  defined  in  that  section  also. 

V  Element  Matrices 

5.1  Introduction 

To  effect  the  discretization  of  the 
element  the  assumed  displacement  functions  are  introduced  into 
the  potential  energy  function  which  in  turn  is  substituted  into 
the  Lagrange  equations  to  yield  element  matrices  with  respect 
to  gridpoint  displacement  degrees  of  freedom.  An  exception  is 
the  element  stress  matrix  which  is  derived  from  strain-displace¬ 
ment  and  stress-strain  relationships. 

5.2  Stiffness  Matrix 

The  energy  contribution  to  the  linear 
elastic  stiffness  is  given  by 

*  =  _lj  [e\  CE]  (e)  dV.  (21) 

2  V 

Substitution  of  the  strain-displacement  relationship. 

Equation  (l8)jln  this  energy  contribution  yields 

$  s  1  f  [Uj  [B]T  [E]  [B]  {U}  dV.  •  (22) 

P  2  J 
V 

Since  matrix  [B]  is  not  a  function  of  the  global  coordinates 
the  Integration  can  be  performed  directly  and 

*P  =  LUJ  CB]T  [E]  [B]  (U) 


U6 


(23) 


The  displacements  {U }  must  he  reordered  to  he  compatible  with 
the  MAGIC  III  ordering  system.  Thus,  the  following  transforma¬ 
tion  is  defined: 


(U>  =  [T]  {U} 


(24) 


where 


U 


z 


1 


* 


The  CT]  matrix  is  defined  hy  Equation  (24A).  Substitution  of 
Equation  (24)  into  Equation  (23)  yields  Equation  (25). 


|p  «  V  LPJ  Ct3t  [B3T  CE]  [B]  Ct]  (u) 


(25) 
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Taking  the  first  variation  of  the  potential  energy  with  respect 
to  the  displacement  degrees  of  freedom  {U >  yields  the  element 
stiffness  matrix  [K]  referenced  to  global  grid  point  displace¬ 
ments.  This  matrix  is  given  below: 

[K]  =  V  [T]T  [B]T  [E]  [B]  [T]  .  (26) 


5-3  Stress  Matrices 

The  element  stress  matrix  follows  as  a  direct 
consequence  of  the  strain-displacement  and  stress-strain  relation¬ 
ships.  Recalling  that 

ia}  =  [E]  [(e)  -  {e}}  (27) 

where  (e)  is  a  column  of  either  mechanical  prestrain  or  thermal 
prestrain  or  both.  Also  recalling  that 

(e)  =  [B]  {U}  =  [B]  [T]  {U>  (28) 

it  follows  that 

{0}  =  [E]  [B]  [T]  (U)  -  [E]  {£}  -  [E]  {e1\  (29) 


ct 

The  vector  {e  },  the  prestrain  due  to  thermal  effects,  can 
be  written  as 


{e°}  =  AT  {a};  {a}  =  |_a  ,  a  ,  a  ,  o,  o,  oj 

a  «y  " 


(30) 


where 


AT  =  T 


ave.  '  V  Tave.  *  T1  +  T2  +  TB  j  TH 

4 


and  Tq  is  a  reference  temperature. 

Equation  (29)  can  now  be  rewritten 
{0}  =  [S]  {UMs}-  {s1}  (31) 
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where: 


CS]  =  [E]  [B]  [T] 

{s}  =  [E]  {e} 

{s>  =  AT  [E]  {a} 

5. 4  Prestrain  Load  Matrix 

The  prestrain  contribution  to  the  potential 
energy  function  is 

=  LeJ  CE]  {e>  dV.  (32) 

Substitution  of  Equations  (18)  and  (24)  into  this  equation 
yields 

^  l°J  ^T]T  [B]T  [E]  lF}  dV 

or  *-  =  v  IAI  CT]T  CB3T  [E]  (F)  =  |uj  .  (33) 

Differentiation  of  Equation  (33)  with  respect  to  the  global 
gridpoint  deformation  yields  the  prestrain  load  vector 

{F-}  =  V  [T]T  [B]T  [E]  (F)  .  (34) 

5.5  Thermal  Load  Matrix 

The  thermal  load  matrix  is  a  special  case  of 
the  prestrain  load  matrix.  Substitution  of  the  thermal  strain, 
Equation  (30),  into  Equation  (34)  yields: 

{Fa>  =  ATV[T]T  [B]T  [E]  {a}  .  (35) 

5.6  Pressure  Load  Matrix 

The  pressure  load  matrix  is  derived  on  the 
basis  of  constant  pressure  on  each  face  of  the  tetrahedron 
element.  Thus  the  total  work,  Wp  }  due  to  the  pressure  loads 
is  the  sum  of  the  work  done  on  each  face. 

Wp  =  W321  +  W431  +  W432  +  W421 
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t 


The  subscripts  denote  a  face  of  the  tetrahedron  (see  Figure 
II-2).  Each  work  term  is  initially  derived  in  a  special  set  of 
local  coordinates  placed  in  a  face  of  the  tetrahedron.  The 
resulting  work  term  is  then  transformed  to  the  global  coordinate 
system.  Thus 


dA 


where  the  negative  sign  accounts  for  direction  of  the 
pressure  and  K  =  s®n  z4* 


(37) 


The  deformation  p  can  be  expressed  in  terms  of  the  assumed 

z 

displacement  functions  and  local  coordinates  x,  y  as 


L V  Mz  » 
h  z2 


where  |yj  =  jx^z^l 


(38) 


Ax  =(-y2Zi|  x  +(x2-x3)  z4^  +  (x3y4“x2y4  +  x4y2  "  x3y2)Z 
+  x3y2z4)  4  |y| 

A2  =(x3z4y  -  x3y4z)  v  |Y| 

A3  =(y2z4  x  “  x2z4y  +  <x2y4  ”  x4y2}  z)* 

a4  =(x3y2z^  M 
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For  purposes  of  integration,  a  triangular  coordinate  system 
is  defined  as  shown  in  Figure  II-3  below: 


The  transformation  from  (x,y)  to  (5,  n)  coordinates  is 
accomplished  by  using  the  following: 


x  =  x2  5  -  (x2-x3)  £n 


y  =  y2(l-n)  5 

dxdy  =  |J(x,y)|d^dn 


|Jtx,y)  |= 


3x 

ay 

ax 

l  =  !-x  V  E  1 

an 

"  an 

as  j  1  x3y2^l 

(39) 


Substitutions  of  Equations  (39)  and  setting  z  =  0  into 
Equations  (38)  yields 


^  *  bv  v  V  V1 


1  -  5 

(l-n)  C 
5n 
o 


(40) 
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Use  of  this  relationship  in  Equation  (37)  and  performing 
the  integration  yields  the  result 


1 
1 
1 
0 

The  work  due  to  pressure  on  face  *<31  is  given  by 
w*<31  *  "  /  p431  py  dA 

A2l3l 


W 


321 


-  ^*tP32i 


l321  |ms 


(Hi) 


(42) 


where  y-  is  the  deformation  parallel  to  the  pressure 

J 

vector.  (See  Figure  II-4  below.) 


y 


y 


As  above. 


the  deformation  y_  can 

v 


be  expressed  as 


=  *-“?■  •  %  •  %  -  ’'yj 


'  1-5 ' 
0 

5n 

(l-n)g 


(43) 


which  is  valid  for  y  =  0.  Thus,  substitution  of  Equation 
(43)  into  Equation  (42)  and  integrating  yields 


>1 


given  by 


The  work  due  to  pressure  on  face  *132  is 


FIGURE  II-5  PRESSURE  LOAD- -FACE  4?2 
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The  yv. ,  i  =  1,  2,  3»  deformations  are  transformed  into 
tetrahedron  local  deformations  using 


{yJ  =  [T1  ]  (y  } 

«7 


(50) 


where 


m 


T 


> 


y 


x 


ll 


e2i>  ^22 *  ®23J  °»  °»  °>  °>  °>  °»  °>  °»  0 

°,  o,  o,  &21,  fc22,  £23,  0,  0,  0,  0,  0,  0 

0,  0,  0,  0,  0,  0,  e21>  e22,  e23,  0,  0,  0 

0,  0,  0,  0,  0,  0,  0,  0,  0  e21,  e22,  e23 

-y2zn 

©22  Zjj (X3  “  ^2^ 

e23  =  y2(X|,  -  x3)  -  y1|(x2  -  x3) 

l%2l  =  C(^21)2  +  (®22)2  +  (®23)2]  1/2 
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The  local  coordinates  are  obtained  fro*  the  transfomabltn 


V 

"V 

'V 

- 

yi 

* 

*i 

- 

Ti 

zi 

1 - 

tor 

L 

1 - 

** 

H 

L. 

( r>0A ) 


1  *  1,2,3^ 


Use  of  Equation  (50)  in  Equation  (49)  gives 


Wl<32  “  “  T  Pji32  A432  CtiJ 


•sT 


0 
1 

i 

1 1 

T 

L1  J 


(51) 


The  work  due  to  pressure  on  face  *121  is  given 


(52) 


where  uy°  is  the  deformation  parallel  to  the  pressure  vector. 
(See  Figure  II-6  below). 


FIGURE  II-6  PRESSURE  LOAD  -  FACE  421 
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The  deformation  jj  j  can  be  written  as 

•J 

wj  -  Lv  °  wy®,  yv  »  j  r  tn 

*"S  v . 

0 

[U-nKj 

which  is  valid  for  y°  »  0*  Substitution  of  Equation  (53) 
into  Equation  (52)  and  integrating  yields: 

W42l  "  *  3  p^31  A^3i  vy2*  ^3*  vy^  T  1  ]  ‘ 


the  u  °  ,  1  *  1,2, 3, 4  deformations  are  transformed  into 
y  i 

tetrahedron  local  deformations  using 


tvyo)  -  [T2]  (?) 


*^ere  "Ifyf*  py0,>  yy°**  yy°*J 

®21*  e22 ,  e23’  ® 
^  0,  0,  0,  ejp  e22*  *23’  °»  01  °*  °*  °*  ° 

^  | e 2 |  ®2l*  ®22 *  e23»  ® 

_°»  0,  0,  0,  0,  0,  D,  0,  0,  e21,  e22>  e23 

e21  *  y2z4,  e*2  --XgZ^,  e*3  -  x2y4~x4y2 


je2!  »  t(e21)2  +  {S22)£  +  (e;3)*r“ 

Use  of  Equation  (55)  in  Equation  (5*0  yields 


2  1/2 


3  pl421  A421  W 
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The  total  work  done  by  a  uniform  normal 
pressure  on  all  sides  of  the  tetrahedron  is  obtained  by  substi¬ 
tuting  Equations  (41),  (46) ,  (51)  and  (56)  into  Equation  (36). 
Thus: 


”,  -  [KtJ •  l^J - 


A.,.*,  cos  0 


-  |  p431  A431 


-  ~  p321  ft321K 


1 

0 

] 

] 

l" 

1 

1 

lOj 


+  ~  p431A431  sin  6 


1/3  P432  A432  t'T^] 

‘  0" 

1 

—  1/3  P421  A421 

”l” 

1 

1 

0 

1 

1 

1 _ 

(57) 


The  local  deformations  ,  U_ 

xi  yi  Zi 


must  first  be  reordered  to  correspond  to  the  MAGIC  III  ordering 
system.  Thus: 

l_b-.fl*  1%J-  L%Jj  ■  [t]T 

In  addition,  these  local  deformations  are  transformed  to 
global  deformations  through  use  of 


(58) 


15]  -  L«J 


(59) 
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VI 


Kinetic  Energy  -  Mass  Matrix 

The  kinetic  energy  for  a  discrete  mechanical 
system  assuming  a  constant  mass  density,  p,  can  be  written: 

*K  *  -f  J  luj  [13  (i>  dv  (61) 

V 

where  £lj  is  an  identity  matrix  and  {y}  are  local  velocities 
of  any  point  in  the  element. 

Lv>J  =  Lv  Uy,  MZJ  •  (62) 

Thus 

*K  =  —  J  [(v^2  +  (yy)2  +  (uz)2]  dv  .  (63) 

v 


The  velocities  y  ,  y  and  y  can  be  expressed  in  terms  of 
a  y  z 

local  grid  point  velocities  through  the  use  of  the  assumed 
displacement  functions.  Thus 


JaL 

6V 


LaJ  =  LA!  j  A2  >  A^ ,  AjjJ 


(6U) 


(65) 
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where : 


Ai  ~  y32zi»  x  +  x23  z2j  y  +  (x3ylj-xi<y3-x2yij+xjjy2 
+  x2y3  -  x3y2)  z  +  6V 

A2  =  "y3z4x  +  x  z4  y  +  (xJly3-x3y1|)  z 
A3  =  y2zi{  x  -  x2z11y  +  (x^-y^ ) 

Aj|  =  (y2x3  -  x2y3)  z 

and  y32  =  y3*-y2  >  x23  =  x2  ”  x3  ' 

Thus  the  products  in  Equations  (63)  are  given  by 

(ux)2  -  Lvx±j  LAj  (vx  } 

(^y)2  =  Li*y±J  {A}  LA1  ^y  }  (66) 

(^z)2  =  Lmz  J  <A)  LAj  {va  >  ,  i  -  1,2, 3, 4  . 


The  kinetic  energy  now  becomes 


(67) 


or 


4k  -  i  L**iJ  E«  <ii> 


where 


CM]  = 


[mO,  0,  O' 


0,  Cm],  0 


t\»m 


0,  0,  [m] 


.[m]  =  pj {A}  [AjdV 


(68) 

(69) 


The  local  grid  point  velocities  in  Equation  (68)  must  be 
reordered  for  use  in  MAGIC  III.  This  is  accomplished 

using  Equation  (58)  in  Equation  (68).  Thus 

*K  "  "!  !xl  [T]T  CM]  [T]  {£}  (70) 


where 

Lfii  -  Lmx  , 

xi 


In  addition,  these  local  grid  point  velocities  must  be 
transformed  to  global  grid  point  velocities.  Equation  (58) 
in  Section  B  of  this  report  is  used.  Thus  the  kinetic  energy  is 

*K  =  -J  lAl  CTg£]T  [T]T  [M3  [T3  [Tg£3  (U)  .  (71) 

Taking  the  first  variation  of  <J>K  with  respect  to  velocities 
and  differentiating  the  result  once  with  respect  to  time 
yields  the  desired  mass  matrix  referenced  to  global  grid 
point  velocities. 
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(72) 


[M]  =  [T gif  [T]T  [M]  [T3  CTg£]  . 

It  now  remains  to  evaluate  the  matrix  Cm] 
of  Equation  (69).  For  purposes  of  integration,  a  tetrahedral 
coordinate  system  will  be  used.  Let  local  coordinates  x,  y, 
z  be  defined  by  the  transformations 

x  =  x4(l-/)  +  x3?f  -  (x3-x2)S'tJ 

y  =  y^d-j')  +  Xjjf?  -  (x3-y2)**if  (73) 

z  =  z4  (l-j') 

dV  =  dxdydz  =  |  J  |  d§d^dj  =  6V  SS2  dsd-^dS*  - 


Using  these  relationships,  the  A  terms  in  Equation  (69)  become 


A1  =  (l~df,  A2  A^  =  >  =  i-j  .  (711) 


Thus  the  integrations  are  performed  simply  and  the  [m] 


matrix  is 


2  111 
12  11 
112  1 
1112 


(75) 


6  3 


D. 


TRIANGULAR  PRISM  ELEMENT 


I .  Introduction 

Three  tetrahedrons  are  assembled  as  shown  in 
Figure  II-7  to  form  a  triangular  prism  element.  When  this 
approach  is  taken,  element  matrices  for  three  tetrahedrons  are 
computed  and  assembled  automatically  within  the  MAGIC  III  System. 
A  considerable  reduction  in  input  is  realized  which  leads  to  a 
corresponding  reduction  in  the  possibility  of  input  error  when 
large  scale  analyses  are  performed.  The  input  required  for  one 
triangular  prism  is  identical  to  that  for  one  tetrahedron  except 
that  six  node  points  define  the  prism  instead  of  four  which 
would  define  the  tetrahedron. 

II .  Element  Static  Matrices 
2 . 1  Stiffness  Matrix 

The  stiffness  matrix  for  each  tetrahedron 
which  makes  up  the  triangular  prism  element  is  computed  in 
accordance  with  Equation  (26)  of  Section  C  of  this  report. 
Recalling  that 

4>p  -  1/2  |0J  [K]  {U}  (1) 

or  for  each  of  the  tetrahedrons: 


Where  Ag  = 


U 


for  example.  The  superscript  I 


refers  to  tetrahedron  number  one.  Also 


4>pZ  =  1/2 

La6>  Al>  a4J 

„II 

K11 

K11 

K12 

K11 

k13 

K11 

Kl4 

A6 

K11 

K22 

K11 

K23 

K11 

K24 

A2 

VII 

K33 

K11 

K3A 

A1 

(Symm. ) 

Ku 

K44 

A4 

Op11  =  1/2 

L^2 »  A5>  a4.1 

“in 

K11 

rrlll 

K12 

JH 

K13 

JII 

Kl4 

'A2  " 

v111 

K22 

nil 

K23 

jh 

K24 

A6 

JII 

K33 

rrlll 

k34 

A5 

(Symm. ) 

Km 

K44 

_V 

The  total  strain  energy  of  the  prism  will  be  the  sum  of  the 
energies  of  the  tetrahedrons.  Thus,  assembly  of  Equations  (2) 
to  (4)  yields: 

4>p  =  1/2  LaPJ  CKP]  {AP}  (5 

P  T 

where  {A  }  =  ^2*  A3*  A45  A5’  A6-* 

and  the  superscript  P  refers  to  prism  quantities, 

p 

[K  ]  is  the  desired  stiffness  matrix  shown  in 
Equation  (6). 


[KP]  = 


+  k”. 

K24  + 

H 

*23J 

k3«' 

K11 

K34 

0 

Kl4  + 

II 

K13 

K22  + 

H 

K22 

k23 

K11  + 
K24  + 

till 

*14 

& 

J-3 

K12  + 

y11 

K12 

+  K111 
+  K1X 

| 

• 

| 

1 

+  k12 

t 

K1  * 

K33  ! 

I 

0 

I 

i 

i 

0 

K13 

(Symmetric) 


II 


k44  +  k44 


lit  III 
K, 


II 


III 


K 


l34 

III 

33 


Kl4  +  K24 


K 


III 

23 


Ku  +  Kn 


+  K 


III 

22 


2.2  Stress  Matrices 
Recalling  that 

{a}  =  [S]  {U}  -  [E]  {e}  -  AT  [E]  (a) 


(7) 


(6) 


Equation  (7)  can  be  used  for  each  tetrahedron  to  give 

to1}  =  [S1]  {U§  -  [E]  {e  I}  -  AT^E]  { a h  (8) 

{a11}  =  [S11]  {U11-  [E]  {i11}-  AT11  [EiHa11}  (9) 

{a111}  =  [SHI]  {U111}  -  [E]  {e111}  -  AT111  [E]{aI3CI}  (10) 
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Neglecting  the  prestrain  and  thermal  strain  for  the  moment. 
Equations  (8)  to  (10)  can  be  rewritten  as: 


[S1I,  s^,  s 3J}  s3i3 


{°oI} 


rq  II  -II  -II  -  II,*-  A 

LS-)  y  y  y  J  A| 


(<,0m)  =  ns,111,  s2m,  s,111,  Sj,111]  ri2 


Assemblying  Equations  (li)  to  (13)  yields  the  desired 
relationship 


rtI  q  I  q  I 

ao  S*»  S2 

II  II  II 

°o  *  S3  S2 


S?  0  0  S1 

0  S  q  g  II 


0  s 


or,  more  compactly; 


III  „  III  S, 


{cP}  =  [SP]  { AP }  . 


The  contribution  to  the  stresses  from 


prestrain  and  thermal  loads  can  be  written  as  follows: 

{Og1}  =  [E]  {e1},  {c^1}  »  [EKe11},  {a£m}=  [E]  {e111}  (16) 


(17) 


or 


{cgP}  =  [EgP]  (eP)  • 

Likewise 

{a^}" AT1  [E]  {a1},  {c*1}  =  AT11  [E]  {a11}, 

(o  m}  =  ATl11  [E]  {a111} 
a 


(18) 


(19) 


(20) 


(21) 
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Combining  equations  (15),  (18),  and  (21)  yields: 


{0P}  =  [SP]  { AP }  -  [E?]  (eP}  -  [EP]  (aP>  •  (22) 


2.3  Pre-Strain  Load 
The  work  energy  is 

<t>£  =  V  [OJ  [T]T  [B]T  [E]  U)  =  [OJ  (P£) 
where 

{F-}=  V  [T]T  [B]T  [E]  {£}  . 
Equation  (23)  can  be  rewritten  as 


where,  for  example: 

LOiJ  =  L»ji,  v  °*iJ  ’ 

Then  for  each  tetrahedron 


(23) 


(2U) 


(25) 
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*<ir> 

e 


La6»  A2’  Al»  A1jJ 


=  a6 ’  A5s  aiJ 


1*1 


2*1 

F3»l 


12,1 


p  1 

1,1 


(ID 


(26) 


F 


2,1 

3,1 


(III) 


(27) 


,P12  *1_ 

Assembly  of  Equations  (25)  to  (27)  yields  the  desired  matrix 

=  l/j  (28) 

where  the  prestrain  load  matrix  is  given  by: 


(F-P 


(I) 

10,1 

+ 

F  UI) 

F7,l 

(I) 

11,1 

+ 

F  (ID 
‘8,1 

(I) 

12,1 

+ 

F  (ID 

*9,1 

(I) 

M 

+ 

77  (H) 

F4,l 

+ 

F  <HD 
Fl,l 

(i) 

5,1 

+ 

77  (II) 

fm 

+ 

f  (HD 

F2,l 

(i) 

6,1 

j 

F  (II> 

F5,l 

+ 

f  (hi) 
F3,l 

} 


7?  <D 
"7,1 


(Continued  on  next  page) 
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P  (ID  (III) 

*  P10,l 

F  (H)  (III) 

11.1  +  FU,1 

p  (H)  (III) 

12-J  +  pi2,i 

Tf  (HI) 

7.1 

p  (HI) 

P8,l 


p  (HI) 

9,1 

P2.‘i>  **2,iW*  F5««> 

F3^  ♦F3.(”>  ♦  >6{?n 


2.H  Thermal  Load  Matriv 


The  prism  thermal  load  matrix  {F  will  be 
of  the  same  form  as  {F?}  in  Equation  (29).  The  force  entries 

4  i.  l.  j  i 


in  this  vector  are  given  by: 


{V  -  VAT  [T]T  [B]T  [E]  (a) 


Equation  (30)  is  evaluated  for  each  tetrahedron. 


i 


2.5  Pressure  Lead  Matrix 

t> 

The  prism  pressure  load  matrix  {F,,  }  will 
be  of  the  same  form  as  {Fg}  in  Equation  (29).  The  force  entries 
in  this  vector  are  given  by : 

{Fp}  -  CTa]T[{CT]T  (Fp)  +  CT1]T  (Fp  )  +  CT23T{$p  >}  *(31) 

1  2 

Equation  (31)  is  evaluated  for  each  tetrahedron. 


Ill .  Kinetic  Energy  and  Mass  Matrix 

The  kinetic  energy  for  each  tetrahedron  is  given 

♦K(I)  -  1/2  lA«,  AiJ  [M(I)]lA«jA*,  As,  AiJT  (32) 

*i11]  -  1/2  U«,Aa,  Ai,  A*  j  [M(II)3  LA,,A*,  Ax,  A,jT  (33) 
*£m)-  1/2  LAj^A.,  As,,  L  J  [M(III)]|Aa,As,  As,  A«,JT  (34) 


The 


mass  matrix  for  the  i 


th 


tetrahedron  can  be 


written  as: 


[M(i)] 


"1.1 

*1,2 

M2,l 

M2,2 

M3.X 

M2,3 

>,1 

M2,4 

Ml,3 

Ml,4 

M2,3 

M2,4 

M3,3 

M3,4 

M3»4 

M4,4 

(35) 


The  total  kinetic  energy  of  the  prism  will  be  the  sum  of  the 
kinetic  energies  of  each  tetrahedron.  Assembly  of  Equations  (32) 
tp  (3*0  yields: 


*P  -  1/2  LAPj  [MP]  (:ap}  (36) 

where 

(AP}T  *  [Ax,  Aa,A»,  A*,  As,  A«J 
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where  Y^,  X2  and  Y2  are  the  global  coordinates  of  the  two 

grid  points  which  define  the  element. 

The  transformation  from  local  deformations  p,  w 
to  global  deformations  U,  V,  W  is  given  by  the  following: 

m  =  [Tgi]  {0}  (2) 

where 


m  »  |vi,  w1#  u2J  w2j 

{U}T  =  luia  V1>  W1,  U 2 »  V2>  W2j 


[T  .]  =  cos  0  sin  80  0 

g  0  0  10 


0  cbs  0  sin8  0 


cos  e  -  — - — -  ,  sin  0 

L 


y2-yi 

L 


III .  Assumed  Displacement  Functions 

The  assumed  displacement  functions  in  the 
local  coordinate  system  are: 


u  =  (ax  +  a2x)  z 

2  3 

w  =  b^  +  b2  x  +  b^  x  +  bjj 


&  •  wx  ■  b2  +  2V  +  3btx2 
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Evaluation  of  Equation  (3)  -a-t  gridpoints  1  and  2  yields: 


- 1 

c 

(-* 

21 

0 

0 

0 

0 

1 

o 

"ar 

W1 

0 

0 

1 

0 

0 

0 

a2 

wx 

X1 

SS 

0 

0 

0 

1 

0 

0 

bl 

u2 

Z2 

X2Z2 

0 

0 

0 

0 

b2 

w2 

0 

0 

1 

X2 

X22 

x  3 
x2 

b3 

c\ 

X 

Se 

0 

0 

0 

1 

2x2 

3x2 

bi|_ 

L.  _ 

— 

_ 

or 

{u}  *  [f]  {A) 


where 

{u}  —  ^1*  wx  *  ^2  *  w2  *  wx2-l 


(4) 


{A}  -  [a^,  q.2»  b^,  b g*  b^,  b(jJ 

“Zl  0  0  0  0 
0  0  10  0 


0 

0 


Thus 


(A>  *  [rr1  {u} 


(5) 
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Evaluation  of  Equatiui.  (3)  gridpoints  1  and  2  yields: 


{u}  =  CP  3  {A}  .  W 

where 


Thus  {A}  =  [r]"1  {u}  (5> 
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Definition  of  the  assumed  displacement  functions 
permits  derivation  of  the  strain-displacement  relationships. 
The  element  strain  components  are: 


or  {e}  =  [B]  {A}  (7) 


(8) 
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where 


Potential  Energy.  Stiffness  Matrix 
The  potential  energy  is 


t/LoJ 


<e)  dV 


{ o }  «  [E]  {e} 


{o)T  -  lcx,  1  ,  axzJ 


[E]  *•  |E  E  0 

L  J  x  xz 


E  E 
zx  z 


0  0 


Use  of  Equations  (5),  (7),  and  (10)  In  Equation  (9)  yields 


*p  -  1/2 


J  LuJ  CP-1;*  £B3T  [E]  [b]  [p”1]  {u}  dV 


*P  *  1/2  [uj  [K]  {u> 


where 


•  /  cr-1; 


3  [B]  [E]  [B]  cr]  dV  . 


Substituting  for  the  [p],  [B]  and  [E]  matrices  from  above 

and  dropping  the  bending  terms  E  .  E_,  E„_  in  the  triple 

X  z  xz 

m 

product  [B]  [E]  [B]  yields  Equation  (13)  after  integration. 
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grid  point . 


The  transformation  matrix  [y]  is  formed  which 


eliminates  the  w  and  w  degrees  of  freedom 
X1  x2 


{u}  *  [y]  {S} 


(15) 


where 


Cy3  «  Cl] 


[K^]"1  [K21] 


.  m  . 

{0}  =  [u^,  w^,  u2,  w2J  and  [I]  is  an  identity  matrix. 


Substitution  of  Equation  (15)  into  Equation  (14)  yields: 

(16) 


♦p  -  -i  |5Jt  :ks]  m 


The  reduced  stiffness  matrix  [KR]  is  given  by 
Equation  (17). 


[Kr] 


CtjT  [%]  [y]  *  AOJ 


where 


-2Z, 


2 

+2Z. 


</*i 

\  l 

-2Z1Z 

"lzT" 


-yi 

V 


\2 


t  G 


xz 


Symmetric 


2’ 


\  Z2  /  . 
2Z 


i!i!  ,VV;2 

L22  I  | 


(17) 


eodiz^  +  38z1z2+  iiz^) 


A  (Zj  +  z2>  CZX2  +  8Zj22  ♦  Z22) 
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,7 

a 


It  is  now  necessary  to  transform  the  local 
deformation  {u}  to  global  deformations.  This  is  accomplished 
by  using  Equation  (2)  in  Equation  (16).  Thus 

*p  -  -i  Luj  CK]  {U>  ( 


where  the  final  desired  stiffness  matrix  for  the  symmetric 
shear  web  is 

tK]  =  CTgt3T  Ckr3  c*  j  . 


V.  STRESS  MATRIX 

In  the  absence  of  prestrain  and  thermal  strain,  the 
stresses  are  given  simply  by  Equation  (10). 

{o}  *  [E]  {e} 

Use  of  Equations  (5)  and  (7)  yields: 


{a}  -  [E]  [B]  IPT1  (u>  • 

The  deformation  vector  {u}  is  reordered  to  be  compatible 
with  Equation  (14)  through  the  transformation 


where 


{u}  =  [Tx]  (u> 


CTj]  - 
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Use  of  Equations  (21), 


(15)  and  (2)  in  Equation  (20)  yields: 


ia)  =  [S]  [T  J  {g} 


(22) 


where  for  x  *  L/2  and  dropping  the  bending  terms  gives: 


[S]  « 


gxz  (6Z12  *  48Z1Z2  +  6Z2}  [~0 
(HZ^  +  38Z.,Z2  +  11Z22)  o 


0  0 
0  0 


G.  TRIANGULAR  RING  ELEMENT  (Asymmetric  Loading) 

I.  Introduction 

The  formulation  of  the  triangular  cross-section 
ring  element  described  herein  is  derived  from,  and  is  mathematically 
consistent  with,  the  formulation  described  in  References  13,  14, 
and  15.  This  ring  element  provides  a  powerful  tool  for  the  analysis 
of  thick-walled  and  solid  axisymmetric  structures  of  finite  length. 
It  may  be  used  to  idealize  any  axisymmetric  structure  taking  into 
account: 

1.  arbitrary  axial  variations  in  geometry, 

2.  axial  variation  in  orientation  of  material 
axes  of  orthotropy, 

3.  radial  and  axial  variations  in  material 
properties, 

4.  any  asymmetric  loading  system  including 
pressure  and  temperature,  and 

5.  degradation  of  material  properties  due 
to  temperature. 


The  discrete  element  technique  was  first 

applied  to  the  analysis  of  axisymmetric  solids  by  Clough  and 

Rashid^ The  formulation  of  the  triangular  cross-section  ring 

(17) 

was  extended  by  Wilson  '  to  include  nonaxisymmetric  as  well 
as  axisymmetric  loads. 

Wilson’s  formulation  for  the  asymmetric  case 

/  n  Q\ 

was  extended  in  Reference  v  '  to  include  orthotropic  material 
properties  with  variable  orientation  axes.  This  extended  develop 
ment  is  presented  here  as  well  as  a  more  precise  means  of  effect- 
ij.ig  the  integration  of  the  strain  energy  over  the  volume  of  the 
ring.  Thermal  and  pressure  load  vectors  and  mass  matrices  are 
also  developed. 

Thus,  the  discrete  element  representation 
presented  consists  of  algebraic  expression  for  the  following 


matrices: 

1. 

Stiffness 

* 

[K] 

2. 

Pressure  Load 

3 

{V 

3. 

Thermal  Load 

3 

<ft> 

4. 

Gravity  Load 

3 

<v 

5. 

Centrifugal  Load 

3 

tcG> 

6. 

Stress 

3 

[S] 

7. 

Mass 

3 

CM] 

The 

matrices  arise  as 

coefficient  matrices  in 

the  Lagrange  equations  for  the  element.  The  appropriate  generalized 
form  of  the  Lagrange  equation  is 

a  ft  +  (  a*2  )  «  o 

3qr  at  3qr 

where 

q^  «  rth  generalized  displacement  coordinate 

■  total  potential  energy 

♦2  »  kinetic  energy 

qr  »  rth  generalized  velocity  coordinate 
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Various  quantities  in  the  following  develop¬ 
ment  will  be  expanded  in  terms  of  Fourier  series .  The  set  of ,  ■ 
unbarred  amplitudes  which  make  up  these  series  are  referred  to 
as  the  A  series  coefficients  and  the  barred  quantities  are 
referred  to  as  the  B  series  coefficients. 

II.  Displacement  Functions  for  the  Triangular  • 

Element 

The  element  generalized  displacements  (see 
Figure  11-10),  can  be  expressed  in  Fourier'  series  form. 

00  - 

u(r,z,8)  =  u  (r,z)  +  I  u. (r,z)  co.s  J8'  +  _?„  u.  (r,z)sinje  (1) 

°  j-1,. J  J*I  ^ 

V(r,z,0)  »  v  (r,z)  +  ?  V. tr,z)  sin  j0  +  I  V\(r,z)cos  je  (2) 
°  J-1  3  ■  j.l  J  . 

88  wo(r,z)  +  2  w,(;r,z)  cos  j0  +  Z  wJrjzJsinje.  (3) 

v  t)  Jar  1 

Linear  displacement  amplitudes  (in  the  r  and  Z  directions) 
are  assumed. 


“j  *  81J  +  62j  r  +  633  z 


VJ  ‘  +  B5J  r  +  66J  z 


"j  *  e7J  +  e8J  r  +  b9J  z 


Note  that  continuity  of  displacement  across  element  boundaries 
is  preserved.  A  transformation  from  generalized  coordinates  to 
grid  point  displacement  coordinates  is  effected  by  writing 


ulj  *  61J  +  62J  V  S3j  Zi 


ViJ  ■  +  B5 1  V  66j  Z1 


wiJ  *  S7J  +  68J  ri  +  s93  Zi 


The  generalized  coordinates,  (BjL> ,  can  be  expressed  (on  the 
harmonic  level)  in  terms  of  grid  point  coordinates  {q^}  as 

{$j}  *  CrBq  3  {<l^}  (8) 

where 


{qj}T  *  ^iy  viy  wij»  U2j *  v2 y  w2j»  u3j»  v3j»  w3jJ  (9) 

{V  *  L*iy  e2 y  B3j»  01»j *  05j*  *6y  B7 y  08j»  e9^  (10) 

Prom  Equation  (7),  with  reference  to  Figure  (TI-10) 


which  is  non-singular. 


(11) 
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Defining  Iqj }  as  follows, 

(qj*>  *  ll'j,  Vj,  Wjj  (12) 

Equations  (4)  thru  (6)  can  be  expressed  in  matrix  form 
as  shown  below 

(qj>  *  u3(r,z)]  {Bj}  •  (13) 

Substituting  Equation  (8)  into  Equation  (13),  an  expression 
relating  the  generalized  element  displacements  to  the  element 
nodal  displacements  (on  the  harmonic  level)  can  be  obtained. 
This  relation  is  given  by  Equation  (14) 

<q*}  =  CX]  {q<J>  (14) 

where 

[X]  =  [$(*»]  [Teq]  •  (15) 


[A]  can  be  expressed  in  explicit  form  as  follows > 


0  0  X2  0  0  x3  0  0 


\1  0 


x2  0 


X1  0 


o  •;  x. 


where 


X1  =  (r2z3"z2r3"(z3"z2)r  +  (r3~r2)Z)/'Al 
X2  *  (zir3"riz3+(z3“zl)r  '■(Vri)  Z)/IA' 

X3  »  (r1z2-z1r2-(z2-z1)r  +  (r^^JZ)/*) A| 
|A|slr  z3  +  rxz2  +  zir3”  Z2V3  “  ri  z3  “  r2zl 
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III.  Potential  Energy 


strain 

The  total  potential  energy  is  derived  as 
energy  and  external  work  contributions. 

the  sum  of 

The  strain  energy  density  is  defined  as 

V  *j"|dej  (a) 

(18) 

where 

- 

(e}T  *  L  err>  ezz»  e00»  exz>  er0»  ez0J 

(19) 

{a}T  *  L°rr»  azz>  a00»  axz»  are>  az0-l  . 

(20) 

Linear  elastic  material  behavior  is  assumed  from 
the  initial  state  of  strain  {ei>  to  the  final  state  of 
stress  {a}  and  strain  {e}, 

(cr(in) }  =  CE(m)]  [{e(m)}  -  ,  (21) 

where  the  superscript  (m)  is  used  to  indicate  that  the 
elastic  modulus  matrix  [E^]  is  evaluated  in  a  coordinate 
system  defined  for  the  material  that  may  be  different  than 
the  r,  z  system  (see  Figure  11-10). 

The  matrix  of  elastic  constants  for  an  orthotropic 
body  with  respect  to  cylindrical  coordinate  axes  is 


Pr(1'V\e)>  W  +  vz9''er>’  Vver  +  V  vea>>  °>  °>  0 
EZ(1-vr8Ver)>  Vv9z  +  vre  °>  °*  0 


1 
A I 


Ee(1-',r2  v*r>>  °>  °>  0 


Symmetric 


AGr3,  0,  0 


AG  „  0 
r0J 
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AG 


z6 


(22) 


i 


where 


A  =  1~vr0v0r'"v0zvz0‘"vzrvrz  “  vr0v0zvzr  “  vrzv0rvz0  .  (23) 


From  symmetry 

Ver  *  E0'lre:  Vzr  *  Ezvrz*  Vez*  Vz6  . 

Poisson’s  ratio,  v^j,  is  defined  as  the  ratio  of  the  strain 

in  the  j  direction  to  the  strain  in  the  i  direction  due  to  a 
stress  in  the  i  direction. 

Equation  (22)  is  more  conveniently  written  in  the  following 
manner : 


(24) 


Substitution  of  the  assumed  constituitive  relations  into  the 
strain  energy  density,  and  the  integration, yields 

u'  =  1/2  [E(m)]  {e(ra)}  -  {e(m)}T  [E(m)]  {e(m)>  (26) 

1  • 

If  the  material  axes  {r^}  are  oriented  at  an  angle  y 
from  the  element  geometric  axes  (see  Figure  11-10),  a  trans¬ 
formation  must  be  introduced 


(e(m))  =  [T£0]  {e} 

(27) 

(o(m)>  »  [Tep]  {a} 

(28) 
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1  ; 


) 

* 

P(m)  P(m) 

E11  E12 

T? 

E13 

0 

0 

0 

\ 

P(m) 

E22 

«(m) 

E23 

0 

0 

0 

t 

! 

> 

p(m) 

E33 

0 

0 

0 

( 

t 

! 

[E(m)]  s 

p(m) 

E44 

0 

0 

(25) 

[ 

f 

I 

4? 

0 

(m) 

b66  J 

i 

■ 

/ 

/ 


cos2y 

sin2y 


2 

sin  y 
cos2y 


0  2slny  cos  y 
0  -2siny  cos  y 


*! 


-siny  cosy  siny  cosy 


2  2 
cos  y-sin  y 


cosy  siny 
-siny  cosy 


Substituting  back  into  Equation  (26)  and  integrating  over 
the  volume  of  the  element,  vre  obtain 

u'-  (1/2  {e}T  [E]  it}  -  {e)T  [E]  {e±})  dV 


where 


[E]  -  [T  ].T  [E(m)]  [T  ]  . 


Equation  (30)  is  the  desired  form  of  the  potential  energy  .■?’< 

The  strains.  Equation  (19)  are  related  to  displacements 
as  follows  in  a  cylindrical  coordinate  system. 

(e}T  -  Lur,  Wz,  U/r  +  VA/r,  U„+W^,  l/r(Uft-V)+VT,,Vz+l/rW J  (32) 


where 


,  Etc. 


Stiffness  Matrix  for  the  Triangular  Element 


In  order  to  effect  the  discretization  of  the  element, 
the  assumed  displacement  functions  are  Introducted  into  the 
potential  energy  function.  Substitution  of  the  total  potential 
energy  function  into  the  Lagrange  equations  yields  the  element 
matrices  with  respect  to  gridpoint  displacements.  Stiffness 
and  mass  matrices,  as  well  as  load  vectors,  are  derived  in 
this  way.  The  element  stress  matrix  Is  derived  from  the  strain- 
displacement  and  stress-strain  relations. 
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The  energy  contribution  of  linear  elastic  stiffness 
is,  in  terms  of  strains, 

*K  -  1/2  {e}T  [E]  {e}  dv  (34) 

In  recognition  of  the  fact  that  the  generalized  dis¬ 
placements  were  described  in  Fourier  series  form,  the  strains 
can  be  described  as  shown  in  Equation  (35) • 

£e>  -  fe  )  +  E  fcjU..}  +  I  fcj  {£.}  (35) 

o  j-i  J  J  J 

For  the  A  series,  j  harmonic  {e^J  is  expressed  as  follows, 

tCJ)T  '  C^j*  ^9,,  c26jJ  <36> 

and  the  Matrix  fc^J  is  a  diagonal  matrix  which  appears  as 
given  in  Equation  (37). 

fcjJ  *  fcos  J  8 ,  cos  j0,  cos  j8,  cos  J 8 ^  sin  J8,  sin  J8J  (37) 
Matrix  | Cj j  is  given  by  Equation  (38). 

[CjJ  ■  fsin  J8,  sin  j 8 ,  sin  J8,  sin  j6,  cos  J8,  cos  J6j  .  (38) 

Expressing  the  strains  (on  the  harmonic  level)  in  terms  of  the 
generalized  coordinates  using  Equations  4,  5,  6.  and  32  yields 

{eJ}  "  fy]  (39) 

where 

m 

{Bj}i  “ 
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and 


CV 


0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

1/r 

1 

Z/r 

J/r 

J 

JZ/r  0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

-J/r 

-J 

-JZ/r 

-1/r  0 

-Z/r  0 

0 

0 

0 

0 

0 

0 

0 

1 

-j/r 

-J 

-JZ/r 

where  for  the  E  series  j  assumes  the  value  of  -j  in  Equation 
(40).  The  differential  volume  is 


(40) 


dV  «  r  d0  dz  dr  . 


(41) 


Substituting  Equations  (35)  and  (4l)  into  Equation  (34), 
and  integrating  with  respects  to  6  yields 

*K  *  2ir  ^  ^  Le#j  r  dz  dr  +  ir£  $^LejKE^ej 

r  z  J*1  T  * 

+  IT  l  ^  ^  L^jJ  EE3  >  rdzdr  . 


}rdzdr 


(42) 


r  z 


It  can  be  seen  that  the  energy  term  represented  by  Equation  (34) 
uncouples  harmonically  (Equation  42)  due  to  the  orthogonality 
conditions  which  exist  mathematically  for  the  triangular  ring. 
The  energy  component  for  the  series,  J*'*1  harmonic  is 


*kj  *  EEfl  rdzdr  ^3) 

r  z 

and  by  substituting  Equation  (39)  into  Equation  (43) 

♦kj  “  ^  UjJ  Cdj3T  EE]  [Dj]  {pj}  r  dzdr  .  (44) 

X*  z 
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Noting  that  the  generalized  coordinates  are  not  variable 
functions  of  r  and  Z,  we  can  write 


rCDj:]T  [E]  [nP  dzdrl  {V  (*5) 

r  z 

where  the  triple  matrix  product  r[Dj]T  CED  [Dj]  is  given  by 
Equation  (46)  on  the  following  page. 

By  inspection  of  the  matrix  in  Equation  (46),  we  see 
that  all  the  integrals  in  Equation  (43)  are  the  type 

6ij  *^ri  Z<J  d2dr  *  w: 

The  integration  is  carried  out  over  the  interior  of 
the  element,  shown  in  Figure  11-10.  The  integration  is 
performed  in  two  parts; 

1)  Between  the  lines  1-2  and  1-3,  i*e.  between 

z  «  k12  r  +  m12  and  z  «  k^  r  +  mi3  from  ri  to  r3* 

2)  Between  the  lines  1-2  and  3-2,  i.e.,  between  z  ■  k12 

r  +  m12  and  z  *  r  +  m^2  from  r3  to  r2* 


where 


K12  * 


Vfl 

r2-rl 


rlZ2"r2Zl 

*2~rl 


Z.-Z, 


r3”rl 


rlZ3~r3Zl 

r3-rl 


Z2~Z3 

r2-r3 


r3Z2~r2Z; 

r2-r3 
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Triple  Matrix  product  r  [D,]a  [E3[D«j 


The  potential  energy  component  for  the  specified  harmonic 
(A  series,  j  harmonic)  is  related  to  the  stiffness  matrix 

for  that  harmonic,  [Kj],  referred  to  generalized  coordinates 
as  follows 

*»  ’  lsi)T  {V  •.  <*> 

[Kj]  is  recognized  as  the  integral  in  Equation  (45).  Its  terms 

are  evaluated  by  substituting  the  appropriate  6^  integrals 

(see  Equation  (47)) for  the  powers  of  r  and  Z  in  Equation  (46) 
as  well  as  the  substitution  of  the  appropriate  harmonic  number  j . 
The  result  is  presented  on  the  following  page  in  Equation  (50). 

Introducing  the  transformation  to  gridpoint  displace¬ 
ments,  Equation  (8)  of  Section  II,  and  taking  the  first  variation 
with  respect  to  the  displacements,  we  obtain  the  element 
stiffness  matrix 

CV  *  [r6%3T  CKp  .  (51) 

Through  a  judicious  choice  of  displacement  functions, 
the  essentially  three-dimensional  character  of  the  ring  changes 
to  one  inherently  two-dimensional  in  nature.  Thus,  an  essentially 
three-dimensional  problem  (asymmetric  loading  on  a  solid  of 
revolution)  can  be  solved  by  undertaking  a  series  of  two-dimen¬ 
sional  applications  of  the  stiffness  matrix  given  by  Equation  (51). 


V.  Load  Vectors  for  the  Triangular  Element 

5.1  Distributed  Load  Vector 

The  external  Work  potential  for  a  system  of  distrib¬ 
uted  loads  (see  Figure  11-10)  acting  on  the  element  face  can 
be  represented  in  the  most  general  form  as  follows: 
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/ 
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Noticing  that  the  generalized  nodal  coordinates  are  not  a 
function  of  ds,  Equation  (6l)  can  be  rewritten  as  follows: 


Vd  “  r  ds  . 

s 


(62) 


Substituting  Equation  (62)  into  the  Lagrange  Equation,  it  can  * 
be  shown  that  generalized  equivalent  nodal  loads  {F^}  can  be 
defined  which  act  on  the  generalized  nodal  coordinates  and  which 
represent  the  mathematical  equivalent  to  the  applied  distributed 
load  system. 


{F  }  can  be  defined  as  follows 

PJ 

{Fp}  =  n^  [X]T  {Pj}r  ds 


(63) 


where 


expressing  the  following  relationships. 


Z  ■  ki2  r  +  M12 


(64) 


where 


K 


12 


Z2  -  Z1 
r2  -  rx 


(65) 


M 


12 


r2Zl  -  rlZ2 
r2  - 


(66) 


and  where 


ds  ■  \/  dr^  +  dz2 


dz 

sin  o( 


(67) 
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it  can  be  shown  that  the  matrix  of  equivalent  nodal  forces 
represented  by  Equation  (63)  can  be  expressed  as  follows 


fe] 


•pie 

,ppj 


K12a 


PrJ{(Al+ClM12)  61  +  (B1  +  cik12>62> 

Pzj{(Al+ClM12)fil  +  <B1  +  C1  k12)6*} 

0 

Prj{(A2  +  C2M12)61  +  (B2  +  C2  K12)62) 


I 

!p30 

n 


PzJ{(A2  +  02M12)«,  +  (B2  +  C2  K12)!2) 


PrJUA3  +  C3M12)5,  ♦  (B3  +  C3  K12)i2) 


PzJ*(A3  +  °3  M12>5‘  +  B3  +  °3  k12) 


The  constants  Ai  and  are  defined  as  follows 


a2 


1 

|A| 

^r2z3”z2r3^ 

B1  * 

1 

|A| 

■’«VZ3> 

1 

1 A  | ' 

(zlr3“rlz3) 

B2  * 

1 

|A| 

-<z3-zl> 

_ 1_ 

-  (r,z0-z,r0) 

B,  » 

1 

•  (z.-z,) 

where 


i*i  *lr2z3 +  viz2  *  2ir3-z2r3  -  ri»z3  -  v-i |  . 

61  and  62  represent  the  following  definite  integrals 

r2 

s“  '  [  r  dr  *  r-2  -  r’2 


(70) 

(71) 


'2  X1 


and 


2  r  3~r  3 

52.  *  \  (jr  «  - - -  • 


S 


(72) 


rl 


A  special  case  is  obtained  when  r2  *  r^.  In  this  instance 
the  formulation  must  be  changed.  For  this  special  case, 
the  equivalent  nodal  load  vector  {F  }  can  be  shown  to  be 


equal  to  Equation  (73). 

f 


».  >  2  < 


Pr,  <Z2  "  zl>  rl 

?z,  (z2  *  zl>  rl 
0 

%  (**-*!>'!  } 
U2-h)rl 

0 
0 
0 
0 


(73) 
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The  load  vectors  represented  by  Equations  (68)  and  (73) 
do  not  account  for  a  distributed  loading  acting  tangentially 

( p 

'  8)  to  the  element  face.  The  loading  system  has  been 
specialized  from  the  original  complete  representation  (Equation 
(53)  to  account  for  varying  distributed  loads  which  act  parallel 
to  the  r  and  Z  axis  (P  and  P  )  respectively.  These  can  be 

combined  to  model  a  varying  pressure  load.  Both  these  loading 
conditions  can  admit  complete  circumferential  asymmetry. 

5.2  Prestrain  and  Thermal  Load  Vectors 

The  prestrain  load  vector  is  constructed  assuming 
uniform  distribution  of  prestrain  across  the  element.  The 
prestrain  contribution  to  the  total  potential  energy  is 

(e}T  [E]  { e ^ }  dV  ,  (7^: 

It  can  be  shown  that  Equation  (52)  when  appropriate  substitu¬ 
tions  are  made  and  an  integratbn  with  respect  to  6  effected, 
takes  the  following  form 

$e  =  2n^  Le0KE3  *ei  >fdzdr  +  11  l*e4jJ[E]{ei>rdzdr 

r  z  0  ^r  z 


+  I  Ej  I  CE]  )rdzdr 


Typically,  for  the  A  series,  j  harmonic  we  have 


■  II  ^  [e^JCE]  it±  Jrdzdr 


i 

! 


Substituting  Equation  (39)  of  Section  IV  into  Equation  (76) 
yields 

=  {VT  CD,]t  rdzdr  [E]  {e±  }  (77) 

j  r  z  J 

Let 

[Dj]  *  n  ^  [d^I1  r  dzdr  (73) 

r  z 


Which  may  be  written  in  terms  of  the 


integrals,  as 


0 

0 

6  0,6 

0 

-J 

^0,0 

6i  »0 

0 

61  0 
y 

0 

-J 

*1,0 

0 

0 

60,1 

*1.0 

-J 

60,1 

0 

0 

0 


[Dj3  - 


0 

0 

0 

0 

0 


0 

0 

0 

G 

6 


J  5 


0,0 


1,0 


J  6 
J  6 
0 
0 
0 


1,0 


0,1 


0 

0 

0 


-6 

0 

-5 


0,0 


0,1 

0 


*1,0  0 


0 

0 

6..  r 

1,0 


0,0 


1,0 

0,lj 
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Transformation  of  Equation  (68)  to  gridpoint  dis 
placement  coordinates  and  substitution  into  the  Lagrange 
equation  yields  the  prestrain  load  vector. 


Where  the  load  components  are 


<Pej}T  *  PfJ.  Py  Plj.  p?>  Plj-  FSjJ  (81> 


and  the  prestrain  components  are 


{eiJ}  *  ^iry  ei®j*  eizj* 


The  thermal  load  vector  is  a  special  case  of  the  pre¬ 
strain  load  vector.  Define  a  matrix  of  thermal  expansion 
coefficients  as 


{a}T  »  |_ar,  aQ,  az>  Oj 


AT  is  the  asymmetric  temperature  rise  above  ambient 
to  which  the  element  is  subjected  and  which  represents  the 
average  of  adjacent  gridpoint  temperatures.  AT  can  be 
expressed  in  Fourier  series  form  as  follows: 

at  »  at  +  v  AT,  cos  J0  +  E  AT,  sin  J6 
o  j*l  i  j*l  J 

i.u  f*  Vi 

The  thermal  load  vector  for  the  J  1  series,  A  harmonic 
appears  as  follows 

{F„  }  *  JI  [IV  ]T  [D.]  CE]  {a}  AT, 


5.3  Gravity  and  Centrifugal  Load  Vectors 

The  external  work  done  by  the  force  of  gravity  on  the 

displacements  can  be  written  as  follows: 

V  =  Cp  GW  dv 
g  J 

dV  =  r  d0  dr  dz 

G  ■  Acceleration  of  Gravity 

p  ■  Mass  Density 

W  ■  Assumed  Displacement  Function 

in  Z  direction 
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Substituting  for  W  into  Equation  (86)  and  integrating 
with  respect  to  8 


Vg  "  211  p°  \  \  Wo  r  drdz 


y  z 


(87) 


Express  Equation  (87)  in  matrix  form  as  follows, 

rdrdz  (88) 


Vg  -  167*,  Sm»  2npG^ 

r  z 


r 

tz 


The  vector  of  forces  on  the  generalized  coordinate  is  in 
terms  of  the  integrals  defined  by  Equation  (47)  of  Section  IV. 
Then 

{Pgo}T  -  2JIpg  [p,  0,  0,  0,  0,  0,  61#,  620  ,  6uJ  (89) 
This  force  is  specifically  a  force  which  is  present  only  in 

i,u 

the  zerotn  or  axisymmetric  harmonic.  The  vector  of  gravity 
forces  on  feridpoint  coordinates  is 


{f  )  -  [r-  ]T  (Fa  ) 


(90) 


The  external  work  done  by  centrifugal  force  due  to  spin  about 
the  axis  of  symmetry  can  be  written  as  follows 

Vs  «$p  u>  r  u  dv  (91) 

where  oi  is  the  spin  rate  and  p  io  the  mass  density,  assumed 
constant  throughout  the  element.  Substituting  for  u  into 
Equation  (91)  and  integrating  with  respect  to  8  gives 


V  *  2IT  pto 
s 


r  dzdr 


(92) 


r  z 
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Expressing  Equation  (92)  in  matrix  form 

vs  *  [$it>  63 ej  2n  pw2  {  r  }r2  dzdr  (93) 

r  z  2 

The  vector  of  forces  as  the  generalized  coordinates  appears  a 

(Fqq)  *  211  poi  ^620,  5 jo  ,  621*  0*  0»  0>  0j  (91*) 

and  the  vector  of  centrigugal  forces  on  gridpoint  coordinates  is 

tPso>  -  (F->  (95) 

Again  {Fgo)  represents  a  force  which  acts  only  on  the 

A. 

zero  or  axisymmetric  harmonic. 

VI .  Stress  Matrices  For  The  Triangular  Element 

The  element  stresses  on  the  harmonic  level  for  the  A 
series,  Jth  harmonic  are  given  by 

{Oj }=  [E]  Uj}  -  [E]  {e^}  (96) 

The  stresses  are  evaluated  at  the  centroid  of  the  cross-section, 
i.e.  at 

rQ  *  1/3  (r^  +  r2  +  r^)  (97) 

ZQ  «  1/3  (z1  +  z2  +  z3> 
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In  Equation  (96),  substitute  for  strains  in  terms  of 
displacements 


{aJ 

}  « 

V  - 

[E]  {e. 
iJ 

} 

(98) 

where , 

from  Equation  (40) 

of  Section 

IV 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

'o 

J 

1 

Zo/ro 

J/ro 

i 

«JZ0/ro 

0 

0 

1 

0 

(99) 

0 

0 

1 

0 

0 

0 

0 

1 

0 

-J/r0 

-j 

-JVro 

-1/r 

0  0 

-Z  /r 
o'  0 

0 

0 

0 

-  0 

0 

0 

0 

0 

1 

-j/ro 

-J 

-JZo/ro 

Equation  (98)  is  used  to  evaluate  elastic  stresses  on  the 
harmonic  level.  The  matrix  {(5}  represents  a  set  of  harmonic 
level  stress  amplitudes.  To  arrive  at  actual  stresses  for  any 
circumferential  position  around  the  element,  the  various  sets 
of  amplitudes  which  arise  during  an  analysis  must  be  recombined 
in  a  set  of  appropriate  Fourier  series.  Thermal  stresses  are 
obtained  by  multiplying  thermal  strains  by  the  matrix  of  elastic 
coefficients.  Equation  (31)  of  Section  III. 

VII .  Mass  Matrix  for  the  Triangular  Element 
The  kinetic  energy  of  the  element  is 

$v  =i^P  (u2  +  V2+w2)  dV  (100) 

1  t  « 

Where  u,  v  and  w  are  the  components  of  radial,  circumferential 
and  axial  velocity.  Substituting  for  u,  v  and  w,  integrating 

with  respect  to  6  and  utilizing  Equations  (8)  from  Section  II, 

4 >v  can  be  cast  in  the  following  form. 
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H.  MODIFIED  QUADRILATERAL  THIN  SHELL  ELEMENT 


I.  Introduction 

The  modified  quadrilateral  thin  shell  element  (Entry  number 
38  in  the  library  of  finite  element  representations  incorporated 
within  the  MAGIC  III  System)  is  described  in  this  section.  This 
finite  element  differs  from  the  present  finite  element  (number  21 
of  the  MAGIC  II  System)  only  in  the  approximation  of  in-plane 
behavior.  No  difference  other  than  the  identification  number  is 
evident  to  the  user. 

This  additional  finite  element  representation,  is  included  in 
the  MAGIC  III  System  for  use  in  the  Idealization  of  membrances  and 
plane-strain  sections  that  require  elongated  finite  element  shapes. 
This  circustance  is  frequently  encountered.  One  important  class 
of  applications  requiring  high  aspect  ratio  finite  elements  is 
the  stress  analysis  of  structural  joints.  A  rule  of  thumb  that 
may  be  applied  to  guide  the  choice  of  element  type  for  such  applica¬ 
tions  is  to  use  the  modified  quadlateral  thin  shell  element  for 
those  finite  elements  whose  aspect  ratio  exceeds  six.  This  guide¬ 
line  derives  from  experience  with  the  IBM  360/65  computer. 

The  approximation  of  in-plane  behavior  embodied  in  the  modified 
quadrilateral  thin  shell  finite  element  differs  from  that  in  the 
original  finite  element  in  two  respects.  Firstly,  the  subdivision 
of  the  finite  element  into  four  triangular  zones  defined  by  the 
diagonals  of  the  quadrilateral  is  avoided  in  generating  the  modified 
finite  element.  This  avoids  the  integrations  over  triangular  zones 
that  were  Judged  to  be  the  principal  constraint  for  accurate  genera¬ 
tion  of  finite  element  number  21  at  high  aspect  ratio.  The  other 
distinguishing  feature  of  the  modified  finite  element  is  that  it 
embodies  a  relatively  simple  discretization  by  direct  interpolation 
of  the  displacement  values  of  the  eight  gridpoints.  The  original 
finite  element  number  21  on  the  other  hand  involves  the  assumption 
of  polynomials  whose  coefficients  must  be  determined  in  terms  of 
the  grldpoint  displacements  by  a  matrix  inversion.  The  accuracy 
of  this  operation  which  is  carried  out  for  each  of  the  four 
triangular  subdivision  deteriorates  with  increasing  aspect  ratio. 
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The  development  and  evaluation  of  the  original  finite 
element  is  presented  in  Pages  113  to  162  of  Reference  1.  The 
development  of  the  modified  finite  element,  number  38,  parallels 
that  of  the  original  finite  element  except  for  the  central  portion 
of  the  representation  of  the  in-plane  behavior.  Therefore,  the 
development  reported  herein  is  confined  to  the  representation  of 
the  in-plane  behavior.  The  interface  of  this  development  with 
that  of  Reference  1  is  clearly  defined  and  a  common  notation  is 
employed.  All  features  of  finite  element  number  21  such  as 
material  orthotropy,  midpoint  node  suppression,  etc.,  are  main¬ 
tained  in  modified  finite  element  (number  38). 

The  implementation  into  the  MAGIC  III  System  leaves  the 
program-analyst  interface  unchanged.  The  user  documentation 
for  finite  element  number  21  applies  to  the  modified  finite 
element  which  is  designated  finite  element  number  38.  The  inter¬ 
face  between  finite  element  library  and  the  surrounding  framework 
of  the  MAGIC  III  System  is  identical  for  finite  element  numbers 
21  and  38.  The  new  calculations  are  confined  entirely  within 
that'  portion  of  the  finite  element  representation  that  generates 
the  basic  in-plane  behavior  representation. 

Numerical  results  are  presented  that  compare  the  original 
and  the  modified  finite  element  representations  at  ordinary  and 
at  high  aspect  ratios.  For  ordinary  aspect  ratios,  the  per¬ 
formance  of  the  modified  finite  element  is  found  to  be  satisfac¬ 
tory  although  generally  less  accurate  than  the  original  finite 
element  which  is  constructed  as  an  assemblage  of  four  subelements. 
However,  for  high  aspect  ratios  the  performance  of  finite  element 
number  38  is  shown  to  be  superior  to  finite  element  number  21. 

This  confirms  the  successful  completion  of  the  effort  to  provide, 
in  the  MAGIC  III  System,  a  quadrilateral  membrane  finite  element 
with  relaxed  constraints  upon  permissible  aspect  ratio. 
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II.  Basic  Relationships 


The  geometry  of  the  quadrilateral  finite  element  is 
illustrated  in  Figure  Il-ll(a).  At  the  branch  point  from  the 
original  sequence  of  calculations  to  the  modified  computation, 
the  following  information  is  known: 


a-  (v  V 

b-  *■ 

c.  CE(gh 

d.  <e(g)> 


coordinates  of  each  of  the  eight 
gridpoints. 

effective  thickness  of  the  membrane. 

material  stiffness  matrix  for  either 
plane  stress  or  plane  strain  as 
appropriate. 

prestress  vector  arising  from  prestrain, 
temperature  load  or  direct  prestress. 


Using  the  foregoing  Information,  the  relations  that  underlie 
the  formulation  of  a  representation  of  the  quadrilateral  membrane 
are  given  below. 

a.  Strain-Displacement  Relation  (Eq.  285,  Ref.  1) 

t_(«)jT  I  auW  „V(E>  au(s>  t  j,v<8) 

|_9,!g  ’  9yg  ’  9yg  9xg 

b.  Stress-Strain  Relation  (Eq.  280,  Ref.  1) 

(o(e>)  -  [|W]((,(«>)  -  (2) 


(1) 


c.  Potential  Energy  Functional  (Eq.  279,  Ref.  1) 

l»(e)j  <3> 

The  construction  of  the  desired  finite  element  representa- 

(e) 

tion  consists  of  the  assumption  of  approximations  for  u  B  and 
v(g)  and  the  substitution  of  these  approximations  into  the  above 
relations.  Then,  integration  of  Equation  (3)  yields  the  basic 
membrane  finite  element  representations,  as: 


i 
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FIGURE  Il-ll(b)  QUADRILATERAL  ELEMENT  IN  TRANSFORMED  SPACE 
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♦m  *  1/2  |>J  CKm<8>5  {V  '  [V]  lF^  W 

<N(g))  =  [Sn(e)3  {6gm>  -  {sn(e)} 
wherein , 

{Sgjjj}  is  the  vector  of  in-plane  gridpoint  displace¬ 
ments  in  the  (x  ,  y  )  coordinate  system  (Eq.  255, 

„  o  o 

Ref.  1) 

[Km^]  is  the  element  membrane  stiffness  matrix  stated 
with  respect  to  the  {6^}  displacement  degrees 
of  freedom. 


{F£  (g)  is  the  element  membrane  prestrain  matrix  stated 
with  respect  to  the  {6^}  displacement  degrees 
of  freedom. 


{N(g)} 


is  the  vector  of  sets  of  membrane  stress  resultants 
aligned  with  the  (x  ,  y  )  coordinate  axes. 

O  o 


is  the  element  stress  matrix  stated  with  respect 
to  the  {5^)  displacement  degrees  of  freedom. 


{sjj^}  is  the  vector  of  sets  of  membrane  prestress 

resultants  aligned  with  the  (xg,  yg)  coordinate 
axes  (Eq.  351,  Ref.  1). 


Equations  (4)  and  (5)  serve  to  define  the  information 
that  the  present  development  must  provide  at  the  point  of 
return  into  the  original  sequence  of  calculations  performed 
in  generating  finite  element  number  21.  Specifically,  the 
matrices  [Km  g^],  {F£^g^}  and  must' be  provided.  The 

vector  {sN(g)}  is  unchanged  by  the  modified  calculations. 


I 


i 

i 

} 


120 


The  present  objective  is  to  develop  explicit  definitions 
for  the  [Km^3»  {Fe^}  and  Once  these  have  been 

obtained,  the  original  sequence  of  calculations  is  reentered 
and  Equations  257*  26l,  262,  263  and  265  of  Reference  1  are 
employed  to  obtain  the  elemenet  stiffness  and  load  matrices 
in  terms  of  the  components  of  displacement  employed  for 
assembly.  This  sequence  of  transformations  can  be  denoted 
symbolically  by: 

* 

Up,,)  >  [r]  <<3>  (5) 

wherein  {q}  is  the  final  set  of  gridpoint  displacement 
degrees  of  freedom.  The  final  form  of  the  finite  element 
I  representation  is  obtained  by  substitution  of  Equation  (6) 

•  into  Equations  (4)  and  (5)  and  adding  to  the  corresponding 

'  representation  of  the  flexural  behavior  in  the  manner  described 

in  Reference.  1. 

I 

\ 

III.  Transformation  of  Coordinates 

1 

;  It  is  clear  from  Equation  (3)  that  the  construction 

I 

!  of  a  finite  element  representation  involves  the  integration 

j  of  functions-  (usually  polynomials)  over  the  interior  region 

of  the  finite  element..  Because  the  performance  of  such 
j  integrations  is  awkward  for  the  quadrilateral  shape  defined 

j  in  the  (x  ,  y  )  coordinates  of  Figure  Il-ll(a)  a  coordinate 

j  0*0 

j  transformation  is  introduced.  Specifically,  the  quadrilateral 

|  element  is  mapped  onto  the  unit  square  of  Figure  II-llb  using 

) 

i  mapping  transformations  defined  by  Reference  20: 


xg(n,y)  *  LH  (n>y)J 

•V 

o<n,u<  1 

(6) 

yg(n,y)  =  LH  (n,y )J 

o<n,y<  1 

(7) 

/ 


wherein: 


<VT  - 

LXel’  Xg2 

vT  • 

[ygl-  :'g2 

S3* 

xg4* 

xg5* 

xg6 

g3* 

yg4. 

gg5» 

yg6 

SV  xg8j 

CO 

S7*  yg8j 

(9) 

{H}  * 


^(l-y)  (+2n  -  2m  -  i) 
ny  (+2n  +2y  -3) 

(l-n)  v  (-2n  +  2y  -  l) 

(l-n)  (l-y)  (~2n  -2y  +  1) 

+4  ny  (l-y) 

+4  ny  (l-n) 

+4  (l-n)  y  (1-fi) 

+4n  (l-n)  (l-y) 

J 


do) 


V 

Given  the  (xg,  y^)  coordinate  values  of  the  eight  grid- 
points,  these  relations  map  the  physical  (Xg,  y  )  space, 
point-by-point,  onto  a  unit  square  in  the  transformed  (n,y) 
space.  Functions  defined  in  the  physical  space  are  expressible 
in  the  transformed  space  as  explicit  functions  of  the  trans¬ 
formed  coordinates,  i.e.. 


f  (x  ,  y  )  «  f  (x  (n,y),  (n,y))  =  f  (n,y)  o<n,y  <  l 

8  8  8  <6 


For  example,  for  the  components  of  displacement  aligned  with 

the  (x  ,  y  )  -  axes: 

©  © 

u(g)  -  u<g)  <„,¥),  v(g)  =  v(K)  0i  n'1'  sl 


Derivatives  of  functions  in  the  (x  .  y„)  coordinates 

8  8 

are  expressible  in  terms  of  derivatives  in  terms  of  the 
transformed  (n,y)  coordinates.  Using  the  chain  rule  of 
differentiation  obtain 
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Tn 
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3y 
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3n 

3x 


8 


3y 


3y 


3n 
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(13) 
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The  inverse  relation  follows  by  direct  calculation,  i.e.. 


r  af  i 

r  syk 

8,k 

r  »f  ^ 

1  X„  I  T 

3y 

3q 

1  g  /  _  1 

7  3f  k 

3x 

8 

4 

3f  \ 

j  3y  | 

{  8 

3y 

3n 

9U 

- 

-J 

^  J 

(14) 


in  which  the  coefficient  matrix  is  denoted  by  [J]  and 
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JQ  =  det  ([J]) 
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(15) 


The  elemental  area  in  the  physical  space  is  related  to  that 
of  the  transformed  space  by: 


dA  =  dXgdy8 


J  d  d 
o  n  u 


(16) 


Equations  (6)  through  (16)  are  sufficient  to  permit 
transformtion  of  the  basic  relations  of  Equations  (1)  through 
(3)  Section  II  to  expression  in  terms  of  the  (n,u)  coordinates. 
The  form  of  the  strain-displacement  relation  becomes, 

fe(g)>  -  [Tu3  UBU}  (17) 


wherein 


{A_n} 

mu 


3u(g)  3u^g) 

~5n  ’  3y 


(18) 
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and: 


[T  3 
L  uJ 


J12» 

0  , 

0 

0 

J21* 

J22 

J22  * 

JH» 

J12 

(19) 


Equation  (17)  is  a  different  mapping  than  that  employed 
in  deriving  finite  element  number  21  but  takes  a  symbolic 
form  identical  to  Equation  299  of  Reference  1. 


As  a  direct  consequence  of  Equation  (17),  the  trans¬ 
formed  stress-displacement  relation  of  Equation  (2),  Section  II 
is  given  by 

{o(g)}  =  [E(g)]  [Tu]  {Amu>  -  (e(g))  (20) 


The  potential  energy  functional  of  Equation  (3),  Section  II 
is  transformed  to  expression  as: 

*  -  cV  (“T  I A  I  [I  t  ]  U  }  dy  (21) 

m  S  \  2  L°muJ  L  mk  mu  -  *-  mu-1  me 

o  o 


wherein 

tv  -  Vo  »„?*  ce(e>J  CV 


(22) 


U 


me} 


t  J  [T  ]T  {e^} 
m  o  u  uJ  * 


(23) 


This  result  is  the  symbolic  equivalent  of  Equation  305 
of  Reference  1  although  the  mapping  employed  is  different. 

The  potential  energy  functional,  as  given  in  Equation  (21), 
is  now  in  a  form  that  readily  admits  integration  over  the 
area  of  the  element  for  the  limits  of  integration  on  ri  and  y 
are  constants. 
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IV.  Discretization 


The  formulation  of  the  finite  element  representation  Is 

carried  forward  by  approximating  the  displacement  functions 
(e)  (*) 

u'B-  and  v'&/  and  integrating  the  potential  energy  over  the 
interior  region  of  the  finite  element.  Polynomials,  defined 
in  the  transformed  space,  are  employed  to  approximate  the  dis¬ 
placement  functions.  The  symbolic  form  of  the  approximations 
is  given  by: 


(g> 

(n,y)  * 

LH 

(n,u)J 

(S(s>) 

(2H) 

(g) 

(n,y)  * 

LH 

<n,u)J 

(25) 

The  vector  of  mode  shapes  (H)  is  the  same  as  that 
employed  to  transform  from  (x  ,  y  )  to  (n,u)  coordinates. 

©  o 

These  mode  shapes  interpolate  the  displacement  functions 
within  the  interior  region  of  the  element  on  the  basis  of  the 
associated  sets  of  gridpoint  displacement  values: 


(26) 


Discretization  of  the  basic  relations  is  accomplished 
in  two  steps.  First,  the  displacement  approximations  are 
employed  to  obtain  {Amu>  of  Equation  (18),  Section  III  as: 

(A  }  *  [D  ]  {6  }  ( 

1  mu  L  nr  gnr 
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wherein : 


3 

LHJ 

>  Loj 

an 

3 

LHJ 

»  Loj 

~3y“ 

L°J 

9 

3 

Sn 

LhJ 

|o| 

9 

3 

|Hj 

— 

3y 

- 

(28) 


(29J 


Now,  using  this  extended  symbolic  notation,  the  basic 
relations  are  discretized.  The  stress-displacement  relation 
of  Equation  (20), Section  III  becomes: 

(„<«*.  [K<g>]  [Tu]  [D„]  (V>-  ISlg)>  <30) 

The  potential  energy  functional  of  Equation  (21),  Section  III 
is  discretized  using  Equation  (27)  to  obtain: 


*»  -  1/2  LVJ  CK<S,]  1 V  -  |>J  {yeCE)}  '3D 


wherein  the  stiffness  [K^]  and  prestrain  load  {F 
matrices  for  the  quadrilateral  membrane  finite  element  are: 


[K9] 


1  1 

Ji 

o  o 


CW  Vu 


(32) 


{P  9>-  1  1  [DlT  (T  }  d  d 


IS 

o  o 


nr  me  n  y 


(33) 
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Two  principal  steps  remain  in  the  development  of  the 
finite  element  representation.  Consideration  must  be  given 
to  the  particularization  of  Equation  (30)  to  specific  points 
within  the  finite  element  and  the  integrations  indicated  in 
Equations  (32)  and  (33)  must  be  carried  out. 

V .  Calculation  of  the  Element  Matrices 

It  is  convenient  to  invoke  numerical  quadrature  to  obtain 
numerical  values  for  the  finite  element  matrices.  All  quantities, 
in  the  integrals  to  be  evaluated  to  obtain  the  element  matrices, 
are  functions  of  the  assumed  mode  shapes  {H}  and  the  gridpoint 
values  {x  }  and  Cy  }.  Thus,  to  obtain  the  values  of  the  inte- 

o  © 

grands,  as  is  required  in  the  numerical  quadrature  calculation, 
it  is  necessary  to  evaluate  the  mode  shapes  {H>  at  the  sample 
points.  Then,  with  these,  numerical  values  can  be  calculated 
for  all  terms  in  the  integrands. 

Gaussian  quadrature  is  employed.  For  the  interval  of 
interest  (of  §fl)  the  set  of  sample  points  (p }  and  weights 
(w)  for  one  dimensional  quadrature  are: 

2- point 

T 

{p}  =  LP.2U32H87,  0.78867483]  (34) 

{oi}  T  *  ip-5  ,  0.5  J  (35) 

3- point 

{p)T  *  JO. 11270165,  0.5  ,  0.887298llJ  (36) 

(w}T  «  [0.27777777,  0.44444444,  0.27777777J  (37) 

These  one-dimensional  sets  of  sample  points  and  weights 

permit  the  construction  of  two-dimensional  sets.  Let  I 

(n,u)>  for  example,  denote  an  integrand  defined  on  the  two- 

dimensional  domain  0  <  n,U  <  1.  Furthermore,  let  the  sample 

points  pr  and  weights  wr  along  the  n  -  coordinate  line  be  R 

in  number.  Similarly,  let  there  be  S  sample  points  p_  and 

weights  wc  along  the  y-  coordinate.  The  Gaussian  product 
s 


127 


formula  for  this  two-dimensional  integration  follows  as: 


The  quadrature  problems  posed  by  Equations  (32)  and  (33) 
Section  IV  involve  integrands  expressed  in  terms  of  {H},  {H,n) 
and  Therefore,  in  preparation  for  quadrature  these 

vectors  are  evaluated  at  the  quadrature  points.  The  collective 
results  are  given  symbolic  definitions  as: 


{H}T 

“  |]AI;q>  1AJi2»  ••• 

Ixl 21  >  L.HJ22* 

• 

•.  LHJ 15 . 

•  »  Dll  25 » 

(39) 

• 

LHJiM»  LhJr2»  ••• 

*» 

{H,n)T 

[LH»nJ  11 »  LH,nJ12, 

*  L?>nJ 21  *  LH*nJ22, 

• 

•  •  • »  LmJ  ^  ^ 

•  •  • »  LH,nJ25 1 

(4o) 

L?»tlIri»  LVl)R2» 

|H,nJR5| 

{H,y}T 

|jLH,y-lii»  Lh>hJi2» 

»  2 1 »  ^>^22* 

1 

L?»hJri»  Lh»hJr2> 

. . . ,  1h,uJ15j 

•  -  • »  25 1 

...»  Lh,u|R5  1 

(41) 

wherein : 


<H,rs 

-  (H  (n,v))| 

<pr,  p8) 

(42) 

=  tH(n,ii))| 

<pr>  P.J 

(43) 

{H,u}rs 

*  Ik  (H(n,u))| 

(pr,  ps) 

(44) 
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The  foregoing  relations  specify  the  quadrature  operation 
completely.  Using  the  evaluated  mode  shapes,  the  element  stiff¬ 
ness  [K^g*]  and  prestrain  load  {Fe^gh  matrices  follow  from 
Equations  (32)  and  (33),  Section  IV  by  direct  calculation,  i.e.. 


“r  “»  CD./  {ImK> 

r*l  s«l 


R  S 

»,<■»  'I  I  “r  “s  CVT  (Ime) 
r«l  s«l 


{H,n)rs,{H,y}rs 


(45> 

(46) 


{H,n)rs,  {H,p)rs 


The  stress-displacement  relation  of  Equation  (30), 

Section  IV  provides  the  means  to  recover  values  for  the 
stresses  at  any  point  within  the  finite  element.  This  rela¬ 
tion  is  particularized  to  a  set  of  five  display  points  similar 
to  that  employed  in  the  original,  number  21,  membrane  finite 
element,  e.g., 

«(e>>!  ■  tv  1  <w  -  1  w> 

wherein : 

[SN(g)]  1  -  [E(g)]  CTU(1,0)]  [Dn(l,0)3  (48) 

(sN}l  -  (e(g)(l,0))  (49) 

The  stress  vectors  at  the  other  points  (n,y)  *  (1,1),  (0,1), 
(0,0)  and  (%,%)  follow  similarly.  The  [Sj^g^]  and  {sN^g^} 
matrices  are  the  matrices  of  Equation  (5),  Section  II  which 
complete  the  specification  of  the  modifications  made  to  the 
original  thin  shell  element  (number  21)  to  obtain  the  modified 
thin  shell  element  (number  38). 
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A  careful  calculation  of  the  gridpoint  loads  that  are 
equivalent  to  a  specified  distribution  of  boundary  loading 
should  be  based  upon  work  equivalence  rather  than  static 
equivalence.  Such  a  calculation  is  not  presently  provided 
within  the  MAGIC  II  System  for  the  membrance  situation  and  must 
be  made  manually.  An  illustrative  calculation  is  included  to 
encourage  the  use  of  work  equivalent  gridpoint  loads. 

Consider  an  element  side  of  Length  L.  For  a  coordinate 
o  s  L  along  the  element  side,  the  assumed  displacement  functions 
employed  in  finite  element  numbers  21  and  30  are  quadratic,  i.e., 

u(s)  =  (s-L)  (s-  j)  uQ 


s  (s-L)  uL^2 


(50) 


Let  the  traction  component  associated  with  this  component 
of  displacement  have  a  specified  distribution,  say  quadratic, 
e.g. , 

p(s)  «  (s-L)  (s-  ^)  pQ 

L 

- —  s  (s-L)  pL/2  (51) 

L 

t  _2_  3  (s-  §>  pL 


The  external  work  of  associated  components  of  boundary 
traction  and  displacement  is  given  by: 


u(s)  p(s)  ds 
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This  energy  functional  is  specialized  to  the  illustrative 
example  by .substitution  from  Equations  (50)  and  (51),  i.e.. 


l. 

V  »  j  |y0,  «L  ,  Uj  r  ^  (a-DCs-  |)"  (s-L)(s-  j>,  -  (s-L)s,  Y?  a  {a-  ^)|p0  ^  dr  (*,3 


-  $  B  <S-L) 
2 
L 

<  ^  B  (S-  Sf 


The  result  of  this  integration  is, 


Luo-  uL/2*  “J 


|o  LU0’  V2>  “J 


8,  i* ,  -2 

32,  4 

-2,  8 


from  which  the  vector  of  gridpoint  loads  is  obtained,  as: 


po  8-  "•  -2  Po 

PL/2  ■  SCI  '•  *  pL/2  ’ 

PL  -2.  1.  8  PL 


This  result  permits  convenient  manual  calculation  of  the 
gridpoint  loads  that  correspond  to  a  quadratic  distribution  of 
boundary  traction  specified  by  its  intensity  at  the  element 
gridpoints. 
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VI .  Convergence 

The  example  chosen  to  illustrate  the  convergence 
characteristics  of  finite  element  number  38  is  the  parabol- 
ically  loaded  membrane  shown  in  Figure  11-12.  This  same 
problem  was  considered  previously  in  evaluation  of  the 
original  finite  element,  number  21. 

The  membrane  is  constructed  of  isotropic  material  and 
the  distributed  loading  is  replaced  by  work  equivalent  grid- 
point  loads  obtained  in  the  manner  outlined  in  Section  V. 

Only  one  quadrant  of  the  membrane  is  considered  explicitly 
in  the  analysis.  This  quadrant  is  idealized  using  square 
finite  elements.  The  idealization  for  the  case  of  four 
finite  elements  is  shown  in  Figure  11-13. 

This  example  prolem  was  analyzed  using  idealizations 
of  1,  4  and  16  finite  elements.  Finite  element  types  21  and 
38  were  considered  as  well  as  a  bi-cubic  element  referred  to 
throughout  as  the  COMEC  finite  element.  Additionally,  a 
solution  obtained  by  an  alternative  method  of  analysis  is 
included  in  the  comparison.  The  displacement  at  the  point  of 
maximum  load,  u  ,  and  the  total  potential  energy  *  are  taken 
to  characterize  the  predicted  behavior. 

The  numerical  results  are  presented  in  Table  II-l.  These 
numerical  results  are  given  graphical  interpretation  in  Figure 
II-14.  It  is  clear  from  Figure  II-14  that  the  maximum  dis¬ 
placement  is  predicted  accurately  by  all  three  types  of  finite 
elements.  Moreover,  the  potential  energy  converges  monotoni- 
cally  for  each  type  of  finite  element.  Specific  displacements 
need  not  converge  monotonically  and  indeed  they  do  not  for  the 
case  illustrated  in  Figure  II-14. 


TABLE  II-l 


PARABOLXCALLY  LOADED  MEMBRANE  CONVERGENCE  RESULTS 


Number  of 
Elements 

Element  Type 

No.  D.O.F. 

Pot.  Energy 

35 

Displace¬ 
ment  u„ 

E  % 
uA 

EXACT 

4 

0.000492 

Magic  Plug  #21 

10 

-0.2138 

0.000492 

0.0 

1 

Magic  Plug  #38 

10 

-0.2147 

0.000496 

0.8 

COMEC 

16 

-0.2162 

0.000489 

0.6 

Magic  Plug  #21 

32 

-0.2155 

0.000492 

0.0 

4 

Magic  Plug  #38 

32 

-0.2167 

0.000492 

0.0 

COMEC 

50 

-0.2169 

0.000493 

0.2 

Magic  Plug  #21 

107 

-0.2167 

0.000492 

0.0 

16 

Magic  Plug  #38 

112 

-0.2169 

0.000492 

0.0 

CO^f' 

162 

-0.2169 

0.000492 

0.0 

j 

|  D.O.F.  =  Degrees  of  Freedom 

i 

t 

! 
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VII.  Shape  Sensitivity 


The  parabolically  loaded  membrane  problem  of  Figure  11-12 
is  employed  to  obtain  an  indication  of  the  sensitivity  of  finite 
element  number  38  to  distortion  of  its  shape  at  ordinary  aspect 
ratios.  The  baseline  idealization  is  comprised  of  four  square 
finite  elements  as  shown  in  Figure  11-13.  Idealizations  of  elements 
of  distorted  shape  are  obtained  by  moving  the  central  gridpoint 
(No.  5)  to  selected  positions  on  the  dashed  circle  shown  in 
Figure  11-13 . 

The  displacement  u  and  the  potential  energy  are  taken. to 
characterize  the  predicted  behavior.  The  results  obtained  using 
finite  element  number  21  are  shown  in  Table  II-2,  together  with 
results  obtained  using  finite  element  number  38  and  the  COMEC 
finite  element.  This  comparison  is  portrayed  graphically  in 
Figure  11-15. 

Observation  of  the  results  of  Table  II-2  and  Figure  11-15 
indicates  that  the  considerable  distortion  imposed  does  not  greatly 
affect  the  accuracy  of  the  behavior  predicted  by  finite  element 
number  38.  It  is  concluded  at  this  point  that  the  new  finite 
element  number  38  may  be  used  in  conjunction  with  the  original 
finite  element  number  21  without  significant  adverse  effects  upon 
the  predicted  behavior.  Indeed,  the  performance  of  the  new  simpler 
finite  element  is  nearly  equivalent  to  that  of  finite  element 
number  21. 

VIII.  Bending  At  High  Aspect  Ratio 

It  is  useful  to  separate  the  evaluation  of  the  performance 
of  finite  element  number  38  at  high  aspect  ratios  into  two  parts. 
First,  bending  is  considered.  Subsequently,  the  type  of  deforma¬ 
tion  which  predominates  in  structural  joints  will  be  examined. 

Consideration  of  bending  at  high  aspect  ratios  is  included 
principally  to  emphasize  the  need  for  caution  in  applications 
where  shear  deformations  are  relied  upon  to  represent  flexural 
behavior.  The  example  problem  chosen  to  illustrate  the  difficulty 
of  coping  with  behavior  of  this  type  is  shown  in  Figure  II-16. 
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TABLE  II-2  PARABOLIC ALLY  LOADED  MEMBRANE  SHAPE  STUDY  RESULTS 


CASE 

ELEMENT  TYPE 

POT.  ENERGY 

S 

DISPLACEMENT 

UQ 

E  % 
u 

EXACT 

- 

0.000492 

- 

1 

Magic  Plug  #21 
Magic  Plug  .#38 
Comec 

-0.2155 

-0.2167 

-0.2169 

0.000492 

0.000492 

0.000493 

0.0 

0.0 

0.2 

2 

Magic  Plug  #21 
Magic  Plug  #38 
Comec 

-0.2164- 

-0.2168 

0.000492 

0.000494 

0.000492 

0.0 

0.4 

0.0 

U 

Magic  Plug  #21 
Magic  Plug  #38 
Comec 

-0.2165 

-0.2166 

0.000491 

0.000492 

0.000489 

m 

— 

4 

Magic  Plug  #21 
Magic  Plug  #38 
Comec 

-0.2166 

-0.2168 

0.000491 

0.000490 

0.000492 

0.2 

0.4 

0.0 

5 

Magic  Plug  #21 
Magic  Plug  #38 
Comec 

-0.2166 

-0.2162 

0.000491 

0.000490 

0.000490 

mm 

6 

Magic  Plug  #21 
Magic  Plug  ^38 
Comec 

-0.2166 

-0.2168 

0.000492 

0.000491 

0.000492 

H 

*  4  FINITE  ELEMENTS 
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,48i_ 


.^8if 


1%  Error 


L  S  MAGIC  PLUG  #21(32  OOP) 
O  =  MAGIC  PLUG  #38(32  OOP) 
□  =  COMEC  (50  DOF) 


CASE 


COMEC 

C 


FIGURE  11-15  PARABOLICALLY  LOADED  MEMBRANE  SHAPE  SENSITIVITY 


The  cantilever  beam  of  Figure  11-16  is  loaded  with  a 
parabolically  distributed  shear  load.  Two  elements,  each  extend¬ 
ing  over  the  entire  depth  are  employed  to  idealize  the  structure. 

A  sequence  of  cases  involving  increasing  aspect  ratios  of  the 
finite  elements  is  obtained  by  holding  the  depth  and  number  of 
finite  elements  constant  while  increasing  the  length  of  the  beam. 

The  displacements,  potential  energy  and  reactions  are 
taken  to  characterize  the  predicted  behavior  of  the  cantilever 
beam.  These  results  are  presented  in  Table  II-3  for  finite 
element  number  38.  Corresponding  results  obtained  from  finite 
element  number  21,  the  COMEC  finite  element  and  beam  theory  are 
also  shown  in  Table  II-3*  Dimensional,  nondimensional  and  error 
values  are  Included. 

Interpretation  of  these  results  is  accomplished  more 
conveniently  using  the  graphical  representation  of  Figure  11-17 . 

At  a  finite  element  aspect  ratio  of  unity,  the  structure  is  not 
a  slender  beam  but  the  finite  element  results  are  in  agreement 
with  each  other  within  a  fraction  of  a  percent. 

At  aspect  ratios  of  two  and  four,  the  finite  element  results 
achieve  reasonable  approximations  of  beam  results  and,  more 
importantly  are  in  satisfactory  agreement  with  each  other  except 
for  the  anomalous  1 error  in  the  reaction  predicted  using 
finite  element  number  38.  The  difficulty  of  representing  bending 
behavior  with  membrane  elements  is  more  apparent  for  the  increased 
aspect  ratio  of  8.0.  While  this  is  not  considered  to  be  especially 
high,  the  system  stiffness  matrix  did  not  admit  accurate  solution 
using  single  precision  arithmetic  on  the  IBM  360/65  computer.  The 
The  reactions  obtained  using  finite  eleirrmt  numbers  21  and  38  are 
grossly  in  error.  Although  not  shown,  for  slightly  higher  aspect 
ratios,  the  reaction  obtained  using  the  COMEC  finite  element  was 
also  grossly  in  error. 
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TABLE  II-3  CANTILEVERED  BEAM,  RESULTS 


DISPLACEMENT  6  =  (U/L3)x  10 


- 


----4k  -> 


1  2  1»  8 

ASPECT  RATIO 

A  =  MAGIC  PLUG  #21  (22  D.O.F.) 

O  =  MAGIC  PLUG  J38  (22  D.O.F.) 

Q  =  COMEC  (37  D.O.F. 

D.O.F.  =  DEGREES  OF  FREEDOM 

FIGURE  11-17  CANTILEVERED  BEAM  BEHAVIOR 
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The  point  of  special  interest  here  is  that  the  source 
of  the  difficulty  does  not  reside  in  the  finite  element  deriva¬ 
tions  themselves.  The  difficulty  is  in  the  conditioning  of  the 
system  stiffness  matrix.  Thus,  the  above  example  emphasizes 
the  inappropriateness  of  this  class  of  finite  elements  for  bending 
applications  but  does  not  constitute  a  meaningful  evaluation 
of  the  relative  performance  of  members  of  this  class  at  high  aspect 
ratios. 

XX.  Tension-Shear  At  High  Aspect  Ratio 

The  results  presented  in'  prior  sections  have  examined 
considerations  that  are  subordinate  to  the  evaluation  of  the 
finite  element  number  38  in  the  present  context.  In  this  section 
of  the  report  the  performance  of  the  modified  finite  element  is 
compared  to  that  of  finite  element  number  21  for  an  idealized 
structural  joint.  Errors  that  arise  in  generating  the  stiffness 
matrix  for  high  aspect  ratio  shapes  of  finite  element  number  21 
have  severely  restricted  attempts  to  analyze  structural  joints 
using  the  IBM  360/65  computer.  The  success  of  the  modification 
of  the  quadrilateral  thin  shell  element  hinges  upon  the  analysis 
of  a  structural  Joint  using  finite  element  shapes  of  substantially 
higher  aspect  ratio  than  Is  possible  with  the  original  finite 
element . 

\ 

The  highly  idealized  structural  joint  employed  in  this 
evaluation  Is  shown  In  Figure  II-18.  Symmetry  permits  explicit 
consideration  of  one  quadrant.  A  tocal  of  four  identical  finite 
elements  arranged  as  shown  in  Figure  11-18  is  used  in  each  case 
considered  in  this  parametric  study.  The  total  load,  uniformly 
distributed  over  the  end,  and  the  length  of  the  Joint  are  held 
constant.  The  parametric  variation  of  the  aspect  ratio  of  the 
finite  elements  is  accomplished  by  varying  the  thickness  of  the 
Joint. 
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The  displacement  uq  on  the  centerline  at  the  load,  the 
potential  energy  and  the  reaction  at  the  line  of  symmetry 
opposing  the  load  are  taken  to  characterize  the  behavior  of  the 
joint.  Of  these,  the  reaction  is  the  most  sensitive  measure. 

The  results  obtained  using  finite  element  number  38  are  compared 
with  those  obtained  using  the  original  finite  element  and  the 
COMEC  finite  element.  Reference  values  are  calculated  consider¬ 
ing  the  Joint  as  a  tensile  bar. 

Two  distinct  series  are  presented  corresponding  to  the 
use  of  isotropic  and  orthotropic  material  properties.  The 
complete  set  of  numerical  results  for  the  isotropic  series  is 
presented  in  Table  II-4.  The  principal  results  are  portrayed 
graphically  in  Figure  11-19 .  It  is  clear  from  Figure  11-19  that 
the  various  predictions  are  in  agreement  at  the  outset.  When 
the  aspect  ratio  is  increased  beyond  8  the  original  finite 
element  representation  leads  to  an  unsatisfactory  error.  On  the 
other  hand,  the  modified  element  representation  performs  satis¬ 
factorily  up  to  a  value  of  64.0.  Thus  the  modified  finite 
element  exhibits  an  improvement  of  a  factor  of  8  over  finite 
element  number  21.  The  relative  accuracy  of  the  COMEC  finite 
element  which  involves  polynomials  of  higher  order  was  unexpected. 

The  same  calculations  were  repeated  for  the  case  of  an 
orthotropic  material.  Table  II-5  contains  the  numerical  results 
and  Figure  11-20  presents  the  corresponding  graphical  Interpre¬ 
tation.  The  original  finite  element  performs  satisfactorily  to 
an  aspect  ratio  of  16,0  while  the  modified  finite  element  is 
apparently  satisfactory  even  beyond  an  aspect  ratio  of  128.0. 

These  results  reinforce  the  factor  of  8  improvement  inferred 
from  the  results  obtained  for  the  isotropic  series. 


TABLE  II-4  ISOTROPIC  LAP  JOINT*  -  RESULTS 
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FIGURE  11-19  ISOTROPIC  LAP-JOINT,  BEHAVIOR 


TABLE  II-5  ORTHOTROPIC  LAP  JOINT  RESULTS 
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FIGURE  11-20  ORTHOTROPIC  LAP  JOINT,  BEHAVIOR 


X.  Summary  and  Conclusions 

The  modification  of  quadrilateral  thin  shell  element 
number  21  was  undertaken  to  relax  the  aspect  ratio  constraint 
on  the  in-plane  portion  of  the  representation.  Attempts  to 
analyze  structural  joints  had  proved  unsuccessful  in  that, large 
residuals  (for  instance,  loss  of  load  throughout  the  structure) 
were  obtained  that  were  attributed  to  the  unavoidable  high 
aspect  ratios  of  the  finite  elements. 

The  development  of  finite  element  number  21  was  examined 
and  the  use  of  triangular  subdivisions  was  Judged  to  be  the 
limiting  factor.  Even  at  modest  aspect  ratios  of  the  quadri¬ 
lateral,  the  ratios  of  the  sides  of  the  triangular  subdivisions 
are  extreme  in  comparison.  Accordingly,  the  principal  modifica¬ 
tion  in  constructing  finite  element  number  38  was  the  elimination 
of  the  use  of  triangular  subdivisions  within  the  finite  element. 

The  modification  to  obtain  finite  element  number  38  is 
presented  in  subsections  I  through  V.  A  simple,  low  order  dis¬ 
placement  approximation  was  chosen  because  experience  has  shown 
that  the  simpler  approximations  are  generally  better  conditioned. 
Additionally,  the  gridpoints  and  grldpoint  degrees  of  freedom  of 
the  final  form  of  the  finite  element  representation  were  stipulated 
at  the  outset  to  be  the  same  as  those  of  finite  element  number  21. 
The  resulting  membrane  representation  of  finite  element  38  is 
equivalent  to  the  quadratic  "serendipity”  isoparametric  finite 
element  representation. 

The  modified  finite  element  representation  is  available  in 
the  MAGIC  III  system  as  finite  element  38.  This  new  finite 
element  representation  maintains  all  features  present  in  finite 
element  21.  The  program-analyst  interface  is  unchanged.  The 
input  data  is  the  same.  The  displays  of  results  have  exactly  the 
same  interpretation  for  the  two  finite  elements. 
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The  results  presented  for  the  membrane  at  different  levels 
of  grid  refinement  establish  that  the  new  membrane  representation 
is  satisfactory  although  somewhat  less  accurate  than  finite  ele¬ 
ment  number  21.  Similarly,  the  results  presented  for  idealiza¬ 
tions  of  the  membrane  using  distorted  finite  element  shapes  show 
that  the  new  element  performs  satisfactorily  at  ordinary  aspect 
ratios . 

The  cantilever  beam  problem  emphasizes  that  this  type  of 
behavior  is  not  predictable  using  full  depth  membranes  (or  shear 
panels)  regardless  of  how  accurately  the  element  matrices  are 
generated.  The  problem  class  of  interest  is  represented  by  the 
idealized  structural  joint  in  which  tension-shear  behavior  is 
dominant . 

The  idealized  isotropic  lap  joint  suggests  an  improvement 
of  a  factor  of  eight  in  the  aspect  ratio  that  can  be  employed 
using  the  new  finite  element  number  38  in  place  of  the  original 
finite  element  number  21.  This  factor  is  substantiated  by  the 
analysis  of  the  same  joint  configuration  constructed  of  ortho¬ 
tropic  materials. 

The  permissible  aspect  ratio  limit  of  finite  element 
number  38  relative  to  finite  element  number  21  is  considered  to 
be  reasonably  well  established  by  the  examples  presented.  How¬ 
ever,  the  permissible  absolute  limit  on  aspect  ratio  depends 
upon  computer  characteristics,  the  size  of  the  problem  and  the 
amount  of  bending  present.  All  results  presented  herein  were 
obtained  using  an  IBM  360/65  computer.  The  problem  sizes  were 
chosen  to  be  small  for  economic  reasons.  Clearly,  the  detri¬ 
mental  effects  of  bending  were  negligible  in  the  illustrative 
lap  joint  examples. 

The  new  quadrilateral  thin  shell  element  number  38  is  re¬ 
commended  for  use  for  elongated  element  shapes  on  the  basis  of 
the  numerical  evaluation  presented  herein.  Its  relative  advan¬ 
tage  is  clear.  Guidance  for  Just  how  large  finite  element  aspect 
ratios  can  be  in  specific  applications  must  evolve  from  usage  in 
practical  design. 


SECTION  III 


INCORPORATION  OP  NEW  COMPUTATIONAL  PROCEDURES 

A.  INTRODUCTION 

Several  new  computational  modules  have  been  Incorporated  Into 
the  MAGIC  III  System  to  support  the  structural  analysis  capability. 

The  first  module  Is  designated  as  ANALIC  (Analysis  In  Core). 
This  module  can  be  used  to  perform  a  complete  linear  elastic  stress 
analysis,  selected  portions  of  a  llnearlly  elastic  analysis  or  as 
a  general  purpose  equation  solver.  Pour  distinct  equation  solvers 
are  available  In  this  module  and  are  described  in  the  following 
subsection.  The  abstraction  instructions  required  for  this  module 
and  detailed  Instructions  for  Its  use  are  delineated  In  Volume  II 
of  this  document  (Reference  7). 

In  addition  to  the  ANALIC  module,  an  additional  out-of-core 
equation  solver  has  been  added  to  MAGIC  III.  A  variable  bandwidth 
solver  utilizing  the  square-root  Cholesky  technique  is  available 
for  the  decomposition  of  symmetric  positive-definite  matrices. 

The  theoretical  details  of  the  method  are  presented  in  a 
subsequent  subsection  while  detailed  instructions  for  its  use  in 
the  MAGIC  III  System  are  given  in  The  User’s  Manual^,  Volume  II 
of  this  report. 
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ANALIC  (ANALYSIS  IK  CORE) 


I,  Introduction 

ANALIC  is  a  MAGIC  III  module  which  can  be  used  to  perform 
a  complete  linear  elastic  stress  analysis  using  in-core  routines. 
This  module  may  also  be  used  to  perform  selected  portions  of  a 
linear  elastic  analysis  or  as  a  general  purpose  equation  solver. 

The  ANALIC  module  is  capable  of  solving  problems  of  approximately 
175  reduced  degrees  of  freedom  with  18,000  words  of  working 
storage.  This  module  features  "dynamic"  storage  which  allows  the 
maximum-sire  problem  to  fit  In  core. 

II.  Equation  Solvers  In  ANALIC 
2.1  Method  of  Bordering 

The  procedure  described  herein  determines  the 
inverse  of  a  symmetric  matrix  by  the  bordering  method.  The  given 
matrix  A  is  regarded  as  the  result  of  bordering  a  matrix  of  order 
(M-l),  the  inverse  of  which  is  assumed  known.  Thus  let 
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Then,  by  seeking  A-1  in  the  same  form,  we  finally  arrive  at 
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where 


The  algorithm  used  is  a  method  of  inverting  a  matrix 
by  successive  borderings.  The  system  loops  on  the  order  of  the 
desired  matrix  inverse  and  computes  the  inverse  of  a  (1  x  1), 
(2x2),  (3x3),  .  .*»  (n  x  n)  in  turn  by  using  the  preceding  com¬ 
puted  Inverse.  Each  step  of  the  process  is  accomplished  on  the 
basis  of  Equation  (2). 

The  following  operations  are  to  be  carried  out  in 

order  to  find  A„  : 

n 

(a)  The  computation  of  the  row  **vnAn~i  elements 

Ynl»  Yn2»  •**»  Yn,n-1 

(b)  The  computation  of  the  number 

n-1 

a  ■  a  +  )  a.  y, 

n  nn  /  ,  in  ’ni 

i-1 


(c) 


The  determination  of  the  elements  alk  of  the 
Inverse  matrix  by  the  relationships 


aik 


i 


(i  <  n-1,  k  <  i) 


(k  <  n-1) 


Storage  for  the  subroutine  used,  consists  of  n  ( — ^ — )  locations 

for  matrix  A  (symmetric  stored  in  the  lower  half  by  rows)  and  one 
column  of  length  n.  The  solution  for  displacements  is  computed 
by  multiplying  the  total  load  column  by  the  computed  inverse. 
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2.2  Gauss  Elimination 

The  subroutine  presented  in  this  section  solves  a 
system  of  simultaneous  linear  equations  with  symmetric  coefficient 
matrix  by  Gauss  elimination.  Consider  the  system  of  simultaneous 
linear  equations 

A  *  X  =  R 

with  symmetric  m  by  m  coefficient  matrix,  the  upper  triangular 
part  of  which  is  stored  by  column  in  m*  (m+l)/2  successive  storage 
locations,  and  an  m  by  n  right-hand  side  matrix  R  stored  by  column 
in  m  *  n  successive  storage  locations .  Solution  is  done  by  Gauss 
Elimination  with  pivoting  in  the  main  diagonal  of  matrix  A.  If 
matrix  R  is  the  identity  matrix,  solution  X  is  the  inverse  of  matrix 
A.  Solution  matrix  X  is  placed  in  positions  of  the  right-hand 
side  matrix  R  and  is  stored  by  column  also.  Thus,  the  computation 
of  the  solution  requires  no  extra  m  by  n  array  of  storage.  Only 
an  auxiliary  storage  array  named  AUX  with  (m-1)  storage  locations 
is  necessary. 

Explicitly,  the  given  system  (1)  is  of  the  formj  + 


+  Note  that  subroutine  GELS  requires  only  the  upper  triangular 
part  of  matrix  Aj  that  is,  the  elements  a..;  a-2;  a22;  a  * 
a33>  •  •  •  5  anm>  a2m»  •••»  amm*  These  elements  are  underlined  in 


“33'  “lm#  2m' 
Equation  (27. 


The  first  step  is  to  search  the  main  diagonal  of  matrix  A  for  the 
element  of  greatest  absolute  value,  say  a^,  and  to  select  it  as 
first  pivot  (p  =  ajj)*  The  reason  ^or  pivoting  only  in  the  main 
diagonal  of  A  is  that  rest-matrices  of  A^  (k  =  1,2, ...,m-l)  must 
remain  symmetrical  during  the  whole  algorithm.  With  a^, 
generate  the  internal  absolute  tolerance  for  testing  usefulness 
of  the  symmetric  algorithm  in  the  following  way: 

tol  =  | adJ J  .6  (3) 

with  a  given  relative,  tolerance  6  . 

Suppose  that  pivot  element  a^  is  equal  to  a^.  if 
it  is  not,  interchange  the  first  rows  of  matrices  A  and  R  with  the 
and  the  first  column  of  matrix  A  with  the  and 'save  column 
interchange  information  by  storing  the  difference  (j-1)  of  pivot 
column  index  j  and  step  counter  k  «  1  [interchanging  column  1 
with  column  j  means  interchanging  variables  x^  with  x^  (1  = 

1,2, ...,n)J  . 

Now  transform  the  elements  of  pivot  rows  in  matrices  A 
and  R  by  division  with  p,  and  the  other  elements  by  adding  -ay  ^  times 
the  new  first  rows  of  these  two  matrices  to  the  other »  rows, 
obtaining ;++ 
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V  -  2, . . . ,m) 

(7) 

-H-  Note  that  transformation  of  pivot  row  in  matrix  A  destroys 
pivot  column,  which  is,  due  to  symmetry,  stored  in  the  same 
location.  As  pivot  column  is  used  unchanged  for  transforma¬ 
tion  of  rest  of  A  and  R,  it  has  to  be  saved  in  auxiliary  array 
AUX  before  transforming  pivot  row. 
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As  column  interchange  information  is  saved  in  the  first  position 
of  the  main  diagonal,  the  result  of  the  first  step  is  the  two 
matrices: 
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It  is  easily  seen  from  equations  (4)  -  (7)  that  the  rest  of  the 
matrix  A^1)  —  that  is,  matrix  A^1^  without  the  first  row  and  first 
column  —  is  symmetric  and  that  actually  only  the  underlined 
elements  must  be  calculated  and  stored.  Therefore,  the  range  of 
index  1  in  formula  (6)  reduces  to  1  =  V,  y  +  1,  . ..,  m* 

This  procedure  is  now  repeated  m~2  times,  starting 
at  each  step  with  the  matrix  of  the  step  before  without  the 
first  k  rows  and  first  k  columns,  and  the  matrix  r(^)  without  the  first 
k  rows.  The  total  result  after  ra-1  steps  is  the  matrices: 
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(Jl-1) 
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Now  work  backward  and  set: 
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(Jo-3)  .. 
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After  each  step  of  back  substitution,  rows  of  solution  matrix 
X  =  r(®)  have  to  be  back- interchanged  according  to  interchange 
information  in  the  corresponding  main  diagonal  element  of  matrix 
A (““l),  in  order  to  get  the  correct  sequence  of  right-hand  side 
column  elements  corresponding  to  the  sequence  of  left-hand  side 
row  elements. 

The  only  case  in  which  the  procedure  described  above  can 
fail  to  give  a  solution  occurs  when  at  any  step  all  elements  in 
the  main  diagonal  of  the  rest-matrix  of  A^)  become  zero,  and  no 
pivot  element  can  be  found.  In  this  case,  the  procedure  is  by¬ 
passed  and  the  error  message  ier  =  -1  is  given.  This  may  — but 
not  necessarily — mean,  that  matrix  A  is  singular.  Possibly 
subroutines  GELG  or  DGELG  (which  are  working  with  complete  pivoting) 
will  be  able  to  find  a  solution  in  cases  where  subroutines  GELS 
or  DGELS  fail.  Actually,  because  of  rounding  errors,  a  further 
check  of  the  absolute  values  of  pivot  elements  is  performed  by 
the  procedure.  If  at  elimination  step  k  this  absolute  value 
becomes  less  than  tol  (see  Equation  3),  it  is  likely  that  there 
was  loss  of  significance  in  the  computation  of  the  diagonal 
elements.  But  as  this  may  not  necessarily  be  the  case,  and  as 
this  test  depends  highly  on  the  choice  of  the  relative  tolerance 
£  ,  the  procedure  gives  only  the  warning  ier  =  k-1,  which  indicates 
that  there  is  a  possible  loss  of  significance  in  the  results 
computed  by  the  algorithm. ++  But  here  it  is  also  possible  that 
subroutines  GELG  or  DGELG  will  give  better  results.  If  there 
is  only  one  equation  to  solve  (m=l),  the  test  on  loss  of  significance 
is  suppressed. 


For  subroutine  GELS,  a  relative  tolerance  between  10 and 
10-L.is  suggested;  and  for  subroutine  DGELS,  between  10 
10 and  10"16. 

For  example,  £  =  10"^  and  warning  ier  =  3  mean  that  there  is  a 
possible  loss  of  about  five  or  more  significant  digits  in  the 
initial  values  of  elimination  step  4. 
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2.3  Cholesky  Triangularizatlon 

Given  an  n  by  n  symmetric  positive  definite  matrix 
A,  compute  an  upper  triangular;  matrix  R  such  that 

A  =  RTR 

The  elements  r^k  of  R  are  computed  using  the 
following  recursive  relationships: 


rlk  “  alk/rll 

rjk  =  (ajk  "  ^  riJrik) 

i=l 


k=l,2,3, . .  .,n 

J~2j 3, ... f n 

k=j, 3+1, . .  .,n 


The  determinant  of  A  is  det(A) 


The  given  matrix  A  is  assumed  stored  columnwise  in 
compressed  form,  that  is  upper  triangular  part  only.  MFSD 
stores  the  solution  R  in  the  same  locations  as  A. 

2 

If  any  calculated  radicand  rkk  (k  =  1,2, ...,n)  is 
not  positive,  further  calculation  is  bypassed,  and  the  error 
parameter  IER  is  set  to  -1.  This  means  that  A  is  not 
positive  definite,  possibly  due  to  roundoff  errors.  IER  Is  also 
set  to  -1  if  the  input  parameter  n  is  less  than  1. 

o 

Let  all  radicands  be  positive  and  let  r^  be  the 
first  radicand  which  is  no  longer  greater  than  the  internal 
tolerance  TOL  =  |EPS  a^J  .  The  subroutine  then  gives  the 
warning  IER  =  k-lj  however,  calculation  is  continued.  The 
warning  indicates  that  there  may  be  loss  of  significance 
at  factorization  step  k  due  to  loss  of  significant  digits  in 

p 

the  calculation  of  rkk. 


Given  a  general  matrix  A  and  a  nonsingular  upper 
triangular  matrix  T,  the  subroutine  MTDS  will  perform  one  of  the 
following  six  operations,  depending  on  the  value  of  an  input 
parameter  IOP: 

IOP=I:-  A  is  replaced  by  T-1A. 

IOP=-l:  A  is  replaced  by  AT"1. 

I0P=2:  A  is  replaced  by  (T-1)TA. 

I0P=-2:  A  is  replaced  by  A(T  )  . 

I0P=3:  A  is  replaced  by  (T  T)  A. 

I0P=-3:  A  is  replaced  by  A(TTT)"1. 

With  the  above  information  available: 

(i)  Calculation  of  X=T_1A  is  done  using  backward 
substitution  to  obtain  X  from  TX=A. 

(ii)  Calculation  of  Y=(T“1)Ta  is  done  using  forward 
substitution  to  obtain  Y  from  TTY=A. 

(iii)  Calculation  of  Z=(TTT)”1A  is  done  by  first 
solving  TTY=A-  and  then  solving  TZ=Y. 

The  remaining  three  operations  are  reducible  to  the 

above  three. 

This  particular  module  may  also  be  used  to  compute 
the  solution  of  a  system  of  equations  BX=A  with  symmetric  positive 
definite  coefficient  matrix  B.  The  first  step  towards  the  solution 
is  the  triangular  factorization  of  B.  The  second  step,  which  may 
be  repeated  for  different  sets  of  righthand  sides  A,  is  the 

rn  l 

calculation  of  (TT)  A.  Another  useful  application  is  the 
computation  of  the  product  ATB~1A  with  symmetric  positive  definite 
B  and  arbitrary  A  in  only  three  steps  and  without  additional 
storage  requirements: 


l6l 


Substituting  into  (2)  we  can  write: 


f 

| 

LUDLUU1  "  f'l  •  K12°2  (5) 


Now  setting 

y  =  dl^ 

we  write 

Llly  =  *1  "  K12°2 

i 

j 

and  solve  for  y  by  forward  substitution.  Finally  we  obtain  the  ! 

unknown  displacements  by  using  backward  substitution  in 
Equation  (7). 

The  stiffness  matrix  is  stored  in  wavefront  format  which  j 
contains  columns  consisting  of  the  first  non-zero  row  to  the  l 

diagonal  element.  The  subroutines  in  ANALIC  operate  on  the  data  I 

in  this  format.  Subroutines  are  called  in  turn  to  convert  the 
symmetric  matrix  to  wavefront  format,  decompose  the  matrix,  perform 
forward  substitution,  and  finally  back  substitution.  j 

| 

i 

j 


(6) 

(7) 
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IV  CONCLUSIONS 


It  is  concluded  that  the  MAGIC  III  System  is  a  logical  and 
consistent  extension  of  the  MAGIC  I  and  II  Systems,  and  that 
the  additional  capabilities  realized  with  MAGIC  III  have  met 
or  exceeded  the  requirements  of  Contract  No.  F33615-71-C-1390. 

The  satisfactory  achievement  of  the  overall  objectives  is 
given  substantiation  by  a  number  of  subsidiary  conclusions. 
Specifically,  it  is  concluded  that: 

(1)  The  addition  of  the  solid  finite  element  represen¬ 
tations  to  the  MAGIC  III  System  provides  enhanced 
capability  to  predict  general  three  dimensional 
states  of  stress  in  structures  of  arbitrary  profile. 

(2)  The  addition  of  the  triangular  cross-section  ring 
finite  element  which  accommodates  asymmetric 
mechanical  and  thermal  loading  on  axisymmetric 
structures  provides  capability  for  the  analysis  of 
thick-walled,  and  solid  axisymmetric  structures  of 
finite  length. 

(3)  The  addition  of  the  modified  quadrilateral  thin  shell 
element  provides  enhanced  capability  for  the 
prediction  of  structural  response  of  membranes  and 
plane-strain  sections  that  require  elongated  finite 
element  shapes. 

(4)  The  addition  of  the  ANALIC  (Analysis  In  Core)  Module 
provides  an  in-core  equation  solution  capability 
designed  for  "moderate-sized"  applications.  Four 
equation  solution  techniques  are  provided. 

(5)  The  out-of-core  variable  bandwidth  equation  solver 
utilizing  the  square  root  Cholesky  technique  has 
been  provided  for  the  decomposition  of  "large  order" 
positive  definite  symmetric  matrices. 
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(6)  The  MAGIC  III  Agendum  Library  has  been  expanded  and 
includes  computational  procedures  for  the  following: 


a.  STATICSASYM 

b.  STATICS 

c.  STATICS C 

d.  STATICS 2 

e.  STABILITY 

f.  STABILITYA 

g.  DYNAMICS 

h.  DYNAMICS F 

i.  DYNAMICS C 

j.  DYNAMICS CF 


(Linear  Elastic  Displacement  and 
Stress  Analysis,  Triangular  Ring  - 
Asymmetric  Loading) 

(Linear  Elastic  Displacement  and 
Stress  Analysis) 

(Linear  Elastic  Displacement  and 
Stress  Analysis  with  Condensation) 
(Linear  Elastic  Displacement  and 
Stress  Analysis  With  Prescribed 
Displacements) 

(Linear  Elastic  Instability  Analysis 
Using  Cholesky  Triangularization) 
(Linear  Elastic  Instability  Analysis 
Using  Matrix  Inversion) 

(Vibration  Frequencies,  Mode  Shapes, 
Generalized  Mass  and  Stiffness  for 
Supported  Structures) 

(Free-Free  Vibration  Frequencies, 

Mode  Shapes,  Generalized  Mass  and 
Generalized  Stiffness  for  Unsupported 
Structures) 

(Vibration  Frequencies,  Mode  Shapes, 
Generalized  Mass  and  Generalized 
Stiffness  with  Condensation  for 
Supported  Structures) 

(Free-Free  Vibration  Frequencies, 

Mode  Shapes,  Generalized  Mass  and 
Generalized  Stiffness  with  Condensa¬ 
tion  for  Unsupported  Structures) 


These  computational  procedures  listed  above  enable 
the  conduct  of  linear  displacement,  stress,  and  - 
stability  analyses  in  the  presence  of  general  prestrain 
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and  thermal  loading  as  well  as  distributed  and 
concentrated  mechanical  loading.  Additionally, 
vibration  analyses  for  free-free  or  supported 
s+ructures  can  be  employed  with  or  without  the 
use  of  condensation  techniques. 

(7)  The  versatile  MAGIC  III  System  finite  element 
library,  which  is  composed  of  sixteen  finite 
elements,  enables  effective  idealization  of  most 
linear  structures. 

(8)  The  stability  analysis  procedure  provided  in  the 
MAGIC  III  System  enables  the  prediction  of 
critical  load  levels  for  general  built-up  shell 
structures. 

(9)  The  preprinted  input  data  forms  facilitate  the  rapid 
and  reliable  specification  of  problem  data  as 
evidenced  by  their  wide  acceptance  with  the 
original  MAGIC  I  and  II  Systems. 

(10)  The  output  provided  by  the  MAGIC  III  System  is 
oriented  to  the  engineering  user,  is  consistent  with 
MAGIC  II,  and  facilitates  clear  and  concise 
interpretation  of  output  parameters. 

(11)  The  computer  program  organization  of  the  MAGIC  III 
Syste  i  is  logical  in  design  and  is  well  suited  to 
generalization. 
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